# Box constraints and weighted sparsity regularization for identifying sources
in elliptic PDEs

Ole Løseth Elvetun and Bjørn Fredrik Nielsen Faculty of Science and
Technology, Norwegian University of Life Sciences, P.O. Box 5003, NO-1432 Ås,
Norway. Email<EMAIL_ADDRESS>of Science and Technology,
Norwegian University of Life Sciences, P.O. Box 5003, NO-1432 Ås, Norway.
Email<EMAIL_ADDRESS>
###### Abstract

We explore the possibility for using boundary data to identify sources in
elliptic PDEs. Even though the associated forward operator has a large null
space, it turns out that box constraints, combined with weighted sparsity
regularization, can enable rather accurate recovery of sources with constant
magnitude/strength. In addition, for sources with varying strength, the
support of the inverse solution will be a subset of the support of the true
source. We present both an analysis of the problem and a series of numerical
experiments. Our work only addresses discretized problems.

The reason for introducing the weighting procedure is that standard
(unweighted) sparsity regularization fails to provide adequate results for the
source identification task considered in this paper. This investigation is
also motivated by applications, e.g., recovering mass distributions from
measurements of gravitational fields and inverse scattering. We develop the
methodology and the analysis in terms of Euclidean spaces, and our results can
therefore be applied to many problems. For example, the results are equally
applicable to models involving the screened Poisson equation as to models
using the Helmholtz equation, with both large and small wave numbers.

Keywords: Inverse source problems, null space, sparsity regularization, box
constraints, PDE-constrained optimization.

##  1 Introduction

We study the task of identifying the source $f$ in an elliptic PDE from
boundary data:

$\min_{(f,u)\in F_{h}\times H^{1}(\Omega),\,0\leq f\leq
s}\left\\{\frac{1}{2}\|u-d\|_{L^{2}(\partial\Omega)}^{2}+\alpha\sum_{i}w_{i}|(f,\phi_{i})_{L^{2}(\Omega)}|\right\\}$
(1)

subject to

$\begin{split}-\Delta u+\epsilon u&=f\quad\mbox{in }\Omega,\\\ \frac{\partial
u}{\partial{\mathbf{n}}}&=0\quad\mbox{on }\partial\Omega,\end{split}$ (2)

where $\Omega\subset\mathbb{R}^{\nu}$ is a Lipschitz domain with boundary
$\partial\Omega$, $\\{\phi_{1},\phi_{2},\ldots,\phi_{n}\\}$ is an
$L^{2}$-orthonormal basis for the finite dimensional function space $F_{h}$,
$\alpha>0$ is a regularization parameter, $d$ represents Dirichlet boundary
data and $\nu=1,2$ or $3$. In this paper we use the same positive weights
$\\{w_{i}\\}_{i=1}^{n}$ in the $\ell^{1}$-regularization term in (1) as in [5,
6]. Their precise form is presented in the next section. The upper bound $s$
in (1) for the strength of the source function $f$ will be either finite or
infinite, i.e., $s\in(0,\infty]$, and we allow the constant $\epsilon$ to be
either positive or negative and of arbitrary fixed size:

$\epsilon=-\kappa^{2}\mbox{ or }\epsilon=\kappa^{2}.$

This means that our results hold for both the screened Poisson equation and
the Helmholtz equation.

Regarding the basis $\\{\phi_{1},\phi_{2},\ldots,\phi_{n}\\}$ for the source
space $F_{h}$, we mention that one can, e.g., for $i=1,2,\ldots,n$, choose
$\phi_{i}$ to be the normalized characteristic function of the $i$’th
cell/subdomain in a mesh for $\Omega$. Provided that the measure of the
overlaps between different cells is zero, this yields an $L^{2}$-orthonormal
basis.

Variants of (1)-(2) occur in a number of applications. For example, in
connection with inverse scattering problems [7, section 1.7]. It has therefore
been studied by many mathematicians, and a number of methods for computing
reliable results have been published.

Often one assumes that $f$ consists of a finite number of point sources or
sources having small compact support, see, e.g., [1, 2, 4, 12, 28]. This leads
to involved mathematical issues, but sometimes optimization procedures or
explicit regularization can be avoided. Also, such ”direct methods” can
potentially be used to identify more general sources [26, 27] when
$\epsilon<0$, i.e., for the Helmholtz equation with multi-frequency data.

In [3, 25] the authors suggest to employ a Kohn-Vogelius fidelity term and a
restricted control domain to identify the source. A similar approach is
considered in [14], but instead of limiting the control domain, a suitable
prior is used.

Assuming that $f$ has the form $f(x)=\rho_{1},\,\rho_{2}$ in
$\Omega_{1},\,\Omega_{2}$, respectively, where $\rho_{1}$ and $\rho_{2}$ are
given constants, the authors of [19, 23] explored the possibility of identify
the subdomains $\Omega_{1}$ and $\Omega_{2}$. These investigations are
motivated by, e.g., the task of determining the mass distribution, i.e. the
core-mantle boundary, of the earth from measurements of the gravitational
field.

If the boundary potential associated with two constant strength sources with
star-shaped supports111Star-shaped with respect to their centers of gravity.
coincide, then these two sources must be identical, see [17, Theorem 4.1.1]
which contains an analysis of this when $\epsilon=0$ (Poisson’s equation). In
[15] a more general case is considered for the Helmholtz equation, and the
author of that paper shows that it is possible to recover the support function
of a polygonal shaped support, see also [16]. Furthermore, during the last
three decades very advanced investigations have been presented to further
clarify the possibilities for computing the support of sources (or scattering
objects) from boundary (or far field) data, see, e.g., [8, 9, 10, 11, 12, 13,
18, 20, 21, 22]. Several of these papers address the important task of
computing the so-called convex source (or scattering) support.

It turns out that the use of box constraints and $\ell^{1}$-regularization in
(1)-(2) enable the recovery of several point/small sources as well as the
identification of rather large composite sources with constant strength,
provided that some fairly mild assumptions hold. Both of these cases can be
handled with the same methodology. We also present results which make it
possible to estimate the constant strength of $f$, as well as the support of
sources with varying strength. Since software for sparsity regularization is
readily available, our scheme can be implemented in a rather straightforward
manner using standard PDE solvers.

Our approach is developed in terms of Euclidean spaces, and it can therefore
be applied to a rather large class of problems. This, e.g., implies that we do
not need to impose limiting restrictions on the domain $\Omega$ and that the
methodology can handle state equations involving conductivity functions or
tensors, i.e., that the term $-\Delta u$ in (2) instead has the form
$-\nabla\cdot(\sigma\nabla u)$.

The mathematical analysis is presented in the next section: Both the basis
pursuit version and the regularized form of the problem are explored. Section
3 contains the numerical experiments, including tests with noisy data, and we
close the paper with some remarks in section 4.

##  2 Analysis

Let

$K_{h}:f\mapsto u|_{\partial\Omega}$

be the forward operator associated with (1)-(2), where $u=u(f)$ denotes the
solution of the state equations (2) for a given $f$. Then we may express our
source identification problem in the form

$\min_{f\in F_{h},\,0\leq f\leq
s}\left\\{\frac{1}{2}\|K_{h}f-d\|_{L^{2}(\partial\Omega)}^{2}+\alpha\sum_{i}w_{i}|(f,\phi_{i})_{L^{2}(\Omega)}|\right\\}.$
(3)

Upon discretization of the boundary value problem (2), we get the
transfer/forward matrix ${\mathsf{A}}\in{\mathbb{R}^{m\times n}}$, and we
obtain the following fully discrete version of (3)

$\min_{0\leq{\mathbf{x}}\leq
s}\left\\{\frac{1}{2}\|{\mathsf{A}}{\mathbf{x}}-{\mathbf{b}}\|_{2}^{2}+\alpha\sum_{i}w_{i}|x_{i}|\right\\}.$
(4)

Here, ${\mathbf{x}}$ represents the Euclidean vector of $f$ relative to the
$L^{2}$-orthonormal basis $\\{\phi_{1},\phi_{2},\ldots,\phi_{n}\\}$ of
$F_{h}$, i.e.,

$f=\sum_{i}x_{i}\phi_{i}.$ (5)

As mentioned above, for the problem (1)-(2) one can, e.g., choose

$\phi_{i}=\|\chi_{\Omega_{i}}\|_{L^{2}(\Omega)}^{-1}\,\chi_{\Omega_{i}},$ (6)

where $\chi_{\Omega_{i}}$ denotes the characteristic function of the $i$’th
cell $\Omega_{i}$ in a mesh for $\Omega$.

Let

${\mathsf{W}}=\mbox{diag}(w_{1},w_{2},\ldots,w_{n})$ (7)

denote the diagonal matrix with diagonal entries $\\{w_{i}\\}$. Then the
associated basis pursuit (zero regularization limit) problem to (4) reads

$\min_{0\leq{\mathbf{x}}\leq
s}\|{\mathsf{W}}{\mathbf{x}}\|_{1}\quad\textnormal{subject
to}\quad{\mathsf{A}}{\mathbf{x}}={\mathbf{b}}.$ (8)

Since only boundary data is available for the recovery of $f$, the operator
$K_{h}$ and the matrix ${\mathsf{A}}$ will have nontrivial null spaces. This
differs significantly from the typical setup for optimal control problems,
which often assumes observations of the state $u$ throughout the entire domain
$\Omega$. Nevertheless, the nature of many important inverse problems is such
that only very limited observation data is available.

Our work is motivated by the inverse problems mentioned above. But for the
sake of generality, we will develop the methodology and analysis in terms of
Euclidean spaces. It should also be mentioned that it remains an open problem
how to establish an infinite dimensional counterpart to the results presented
in this paper.

###  2.1 Diagonal weight/regularization matrix ${\mathsf{W}}$

Assume that ${\mathsf{A}}\in{\mathbb{R}^{m\times n}}$ has a nontrivial null
space ${\mathcal{N}({\mathsf{A}})}$ and consider the equation

${\mathsf{A}}{\mathbf{x}}={\mathbf{b}}.$ (9)

We can multiply (9) with the pseudo-inverse ${\mathsf{A}^{\dagger}}$ of
${\mathsf{A}}$ to obtain

${\mathsf{A}^{\dagger}}{\mathsf{A}}{\mathbf{x}}={\mathsf{A}^{\dagger}}{\mathbf{b}},$
(10)

where we recall that ${\mathsf{A}^{\dagger}}{\mathsf{A}}$ is the orthogonal
projection

$\mathsf{P}={\mathsf{A}^{\dagger}}{\mathsf{A}}:{\mathbb{R}^{n}}\rightarrow{\mathcal{N}({\mathsf{A}})^{\bot}}.$
(11)

The authors of [5] introduced the diagonal weight/regularization matrix
${\mathsf{W}}\in{\mathbb{R}^{n\times n}}$ defined by222In section 3 in [5],
${\mathsf{W}}$ is defined in terms of two abstract finite dimensional inner
product space $X$ and $Y$ and a linear operator ${\mathsf{K_{h}}}:X\rightarrow
Y$. In this paper we use the ”Euclidean version” of this.

${\mathsf{W}}{\mathbf{e}}_{i}:=\underbrace{\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}}_{:=w_{i}}{\mathbf{e}}_{i},\quad
i=1,2,\ldots,n,$ (12)

where ${\mathbf{e}}_{1},{\mathbf{e}}_{2},\ldots,{\mathbf{e}}_{n}$ are the
standard unit basis vectors. Throughout this paper we assume that none of the
standard unit basis vectors belong to the null space
${\mathcal{N}({\mathsf{A}})}$ of ${\mathsf{A}}$, which implies that

$w_{i}=\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}\neq 0,\quad i=1,2,\ldots,n,$ (13)

and ${\mathsf{W}}$ becomes invertible since it is diagonal (7).

Some of the mathematical properties of ${\mathsf{W}}$ are explored in [5, 6].
In order to motivate the use of this type of weighting procedure, we mention
the following: Let ${\mathbf{x}}^{\dagger}$ be the minimum
$\|\cdot\|_{2}$-norm solution of

${\mathsf{A}}{\mathbf{x}}={\mathsf{A}}{\mathbf{e}}_{j}.$ (14)

Then, see Theorem 2.1 in [6] and Theorem 4.2 in [5],

$\displaystyle j$
$\displaystyle=\operatorname*{arg\,max}_{i\in\\{1,2,...,n\\}}|[{\mathsf{W}}^{-1}{\mathbf{x}}^{\dagger}]_{i}|$
$\displaystyle=\operatorname*{arg\,max}_{i\in\\{1,2,...,n\\}}|[{\mathsf{W}}^{-1}{\mathsf{A}}^{\dagger}{\mathsf{A}}{\mathbf{e}}_{j}]_{i}|,$
(15)

where $[\mathbf{v}]_{i}$ represents the $i$’th component of the Euclidean
vector $\mathbf{v}$. This means that ${\mathsf{W}}^{-1}{\mathbf{x}}^{\dagger}$
can be used to recover the correct index $j$ from the ”data”
${\mathbf{b}}^{\dagger}={\mathsf{A}}{\mathbf{e}}_{j}$ in (14).

#### Remarks

In order to obtain a ”sound” regularization operator ${\mathsf{W}}$, as
defined in (12), one should ensure that
$\min_{i}\\{\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}\\}$ is not too small. This
follows from the observation that, if
$\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}\approx 0$, then one tries to recover the
contribution associated with a basis vector ${\mathbf{e}}_{i}$ which is almost
in the null space of the matrix ${\mathsf{A}}$. Thus, in order to avoid to get
”very close” to the null space, one would typically use a rather moderate
number $n$ of basis functions $\phi_{1},\phi_{2},\ldots,\phi_{n}$ to
discretize the source $f$ in (1)-(2), see (5) and (6). This is our motivation
for employing a finite dimensional space for $f$ in (1)-(2).

The computation of
$\mathsf{P}{\mathbf{e}}_{i}={\mathsf{A}}^{\dagger}{\mathsf{A}}{\mathbf{e}}_{i}$
involves the pseudoinverse/Moore-Penrose inverse ${\mathsf{A}}^{\dagger}$ of
${\mathsf{A}}$. If ${\mathsf{A}}$ has small nonzero singular values, which
will typically be the case for a discretized ill posed problem, then it is not
advisable to apply ${\mathsf{A}}^{\dagger}$ in computations, and we must
”replace” ${\mathsf{A}}^{\dagger}$ by a more well-behaved matrix. This can,
e.g., be accomplished by employing truncated SVD or by applying Tikhonov
regularization, with regularization parameter $\gamma$,

$({\mathsf{A}}^{T}{\mathsf{A}}+\gamma{\mathsf{I}}){\mathbf{x}}_{\gamma}={\mathsf{A}}^{T}{\mathsf{A}}{\mathbf{e}}_{i}.$

Then, according to standard theory, $\lim_{\gamma\rightarrow
0^{+}}{\mathbf{x}}_{\gamma}={\mathsf{A}}^{\dagger}{\mathsf{A}}{\mathbf{e}}_{i}=\mathsf{P}{\mathbf{e}}_{i}$.
In our numerical experiments (section 3) we used truncated SVD to compute an
approximation of the entries in the diagonal regularization matrix
${\mathsf{W}}$.

###  2.2 Basis pursuit

We will now investigate whether one can use333If ${\mathsf{A}}$ has very small
nonzero singular values, then it is in practise not possible to solve (8), and
one must instead consider a regularized version of the problem. The latter is
discussed in the next subsection. (8) to approximately, or partially, recover
a true source

${\mathbf{x}}^{*}=\sum_{\mathcal{J}}x_{j}^{*}{\mathbf{e}}_{j},\quad{\mathbf{x}}^{*}\geq
0,$

from the synthetic data
${\mathbf{b}}={\mathbf{b}}^{\dagger}={\mathsf{A}}{\mathbf{x}}^{*}$, where
$\mathcal{J}$ is the support of ${\mathbf{x}}^{*}$.

Our first result shows that requiring positivity may enable the identification
of the support of the sources.

######  Theorem 2.1.

Let $\mathcal{J}=\mathrm{supp}({\mathbf{x}}^{*})$ and assume that there exists
a vector ${\mathbf{c}}$ such that

$\displaystyle\frac{\mathsf{P}{\mathbf{e}}_{j}}{\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}}\cdot{\mathbf{c}}$
$\displaystyle=1\quad\forall j\in\mathcal{J},$ (16)
$\displaystyle\frac{\mathsf{P}{\mathbf{e}}_{i}}{\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}}\cdot{\mathbf{c}}$
$\displaystyle<1\quad\forall i\in\mathcal{J}^{c}.$ (17)

Then any solution ${\mathbf{y}}$ of

$\min_{0\leq{\mathbf{x}}\leq
s}\|{\mathsf{W}}{\mathbf{x}}\|_{1}\quad\textnormal{subject
to}\quad{\mathsf{A}}{\mathbf{x}}={\mathsf{A}}\overbrace{\sum_{\mathcal{J}}x_{j}^{*}{\mathbf{e}}_{j}}^{={\mathbf{x}}^{*}}$
(18)

satisfies

$\mathrm{supp}({\mathbf{y}})\subseteq\mathrm{supp}({\mathbf{x}}^{*}),$

provided that ${\mathbf{x}}^{*}\geq 0$ and that
$s\in[\|{\mathbf{x}}^{*}\|_{\infty},\infty]$. Also, ${\mathbf{x}}^{*}$ is a
solution of (18).

#### Remark

Our geometric intuition for assumptions (16) and (17) is: There exists a
direction ${\mathbf{c}}$ such that
$\mathsf{P}{\mathbf{e}}_{j}/\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}$, $\forall
j\in\mathcal{J}$, and
$\mathsf{P}{\mathbf{e}}_{i}/\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}$, $\forall
i\in\mathcal{J}^{c}$, are ”equally important” and ”less important”,
respectively, in this direction. We will in subsection 2.3 explain how (16)
and (17) can be derived from the classical KKT-conditions.

###### Proof of Theorem 2.1.

Let

${\mathbf{y}}=\sum_{\mathcal{J}^{c}\cup\mathcal{J}}y_{i}{\mathbf{e}}_{i},\quad
y_{i}\geq 0,$

denote a solution of (18). Then

${\mathsf{A}}{\mathbf{y}}={\mathsf{A}}\sum_{\mathcal{J}}x_{j}^{*}{\mathbf{e}}_{j}$

or, because $\mathsf{P}={\mathsf{A}}^{\dagger}{\mathsf{A}}$,

$\mathsf{P}{\mathbf{y}}=\mathsf{P}\sum_{\mathcal{J}}x_{j}^{*}{\mathbf{e}}_{j}.$

Taking the inner product with ${\mathbf{c}}$ yields

$\sum_{\mathcal{J}^{c}\cup\mathcal{J}}y_{i}\mathsf{P}{\mathbf{e}}_{i}\cdot{\mathbf{c}}=\sum_{\mathcal{J}}x_{j}^{*}\mathsf{P}{\mathbf{e}}_{j}\cdot{\mathbf{c}}$

which can be written in the form, using (16),

$\sum_{\mathcal{J}^{c}\cup\mathcal{J}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}\frac{\mathsf{P}{\mathbf{e}}_{i}}{\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}}\cdot{\mathbf{c}}=\sum_{\mathcal{J}}x_{j}^{*}\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}.$
(19)

Assume now that there is $i^{\prime}\in\mathcal{J}^{c}$ such that
$0<y_{i^{\prime}}\leq s$. Then, combining (19) with (16), (17) and the
inequality constraint $y_{i}\geq 0$ lead to the conclusion that

$\sum_{\mathcal{J}^{c}\cup\mathcal{J}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}>\sum_{\mathcal{J}}x_{j}^{*}\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}.$
(20)

That is, recall the definition (12) of ${\mathsf{W}}$,

$\displaystyle\|{\mathsf{W}}{\mathbf{y}}\|_{1}$
$\displaystyle=\|\sum_{\mathcal{J}^{c}\cup\mathcal{J}}y_{i}{\mathsf{W}}{\mathbf{e}}_{i}\|_{1}$
$\displaystyle=\sum_{\mathcal{J}^{c}\cup\mathcal{J}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}$
$\displaystyle>\sum_{\mathcal{J}}x_{j}^{*}\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}$
(21) $\displaystyle=\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1},$

where we have used the assumption that ${\mathbf{x}}^{*}\geq 0$ in the last
equality. Consequently, since ${\mathbf{x}}^{*}$ satisfies the constraints in
(18), there can not be any $i^{\prime}\in\mathcal{J}^{c}$ such that
$y_{i^{\prime}}>0$. We conclude that
$\mathrm{supp}({\mathbf{y}})\subseteq\mathcal{J}$.

When $\mathrm{supp}({\mathbf{y}})\subseteq\mathcal{J}$, we must replace $>$
with $\geq$ in (20) and (21). It follows that

$\|{\mathsf{W}}{\mathbf{y}}\|_{1}\geq\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}$

for any solution ${\mathbf{y}}$ of (18), and we conclude that
${\mathbf{x}}^{*}$ solves (18). ∎

As an immediate consequence we obtain the following result, where we remark
that: Gradually decreasing $s$ in (18) may enforce uniqueness because the size
of
$\left\\{{\mathbf{z}}\,|\,\mathrm{supp}({\mathbf{z}})\subseteq\mathcal{J}\mbox{
and }0\leq{\mathbf{z}}\leq s\right\\}$ will decrease as $s$ decreases.

######  Corollary 2.1.1.

Assume that ${\mathbf{x}}^{*}=\sum_{\mathcal{J}}x_{j}^{*}{\mathbf{e}}_{j}$,
${\mathbf{x}}^{*}\geq 0$, is such that a vector ${\mathbf{c}}$ satisfying (16)
and (17) exists. If ${\mathsf{A}}$ is injective on
$\left\\{{\mathbf{z}}\,|\,\mathrm{supp}({\mathbf{z}})\subseteq\mathcal{J}\mbox{
and }0\leq{\mathbf{z}}\leq s\right\\}$, where
$s\in\left[\|{\mathbf{x}}^{*}\|_{\infty},\infty\right]$, then
${\mathbf{x}}^{*}$ is the unique solution of (18).

Recall the definition of the ”$\ell^{0}$-norm”

$\|{\mathbf{v}}\|_{0}=\\#\\{i\,|\,v_{i}\neq 0\\}.$

Roughly speaking, the next result explains that, if a minimum
”$\ell^{0}$-norm” solution with support $\mathcal{J}$ obeying (16) and (17)
exists, then all solutions of the basis pursuit problem must have the same
support – i.e., the support is unique.

######  Corollary 2.1.2.

Assume that

${\mathbf{b}}\in\\{{\mathsf{A}}{\mathbf{z}}\,|\,{\mathbf{z}}\geq 0\\},$

and that there exist

$\widetilde{{\mathbf{x}}}^{*}\in\operatorname*{arg\,min}_{{\mathbf{x}}\geq
0}\|{\mathbf{x}}\|_{0}\quad\textnormal{subject
to}\quad{\mathsf{A}}{\mathbf{x}}={\mathbf{b}}$ (22)

and a vector ${\mathbf{c}}$ such that (16) and (17) hold with
$\mathcal{J}=\textnormal{supp}(\widetilde{{\mathbf{x}}}^{*})$. Then any
solution ${\mathbf{y}}$ of

$\min_{0\leq{\mathbf{x}}\leq
s}\|{\mathsf{W}}{\mathbf{x}}\|_{1}\quad\textnormal{subject
to}\quad{\mathsf{A}}{\mathbf{x}}={\mathbf{b}}$ (23)

also solves the problem

$\min_{0\leq{\mathbf{x}}\leq s}\|{\mathbf{x}}\|_{0}\quad\textnormal{subject
to}\quad{\mathsf{A}}{\mathbf{x}}={\mathbf{b}}$

and

$\mathrm{supp}({\mathbf{y}})=\mathrm{supp}(\widetilde{{\mathbf{x}}}^{*}).$

This holds for any
$s\in\left[\|\widetilde{{\mathbf{x}}}^{*}\|_{\infty},\infty\right]$.

###### Proof.

Since ${\mathsf{A}}\widetilde{{\mathbf{x}}}^{*}={\mathbf{b}}$, any solution
${\mathbf{y}}$ of (23) must also be solution of (18) with
${\mathbf{x}}^{*}=\widetilde{{\mathbf{x}}}^{*}$. Therefore, according to
Theorem 2.1, ${\mathbf{y}}$ satisfies

$\mathrm{supp}({\mathbf{y}})\subseteq\mathrm{supp}(\widetilde{{\mathbf{x}}}^{*}).$
(24)

Consequently,

$\|{\mathbf{y}}\|_{0}\leq\|\widetilde{{\mathbf{x}}}^{*}\|_{0},$

but from the definition (22) of $\widetilde{{\mathbf{x}}}^{*}$

$\|{\mathbf{y}}\|_{0}\geq\|\widetilde{{\mathbf{x}}}^{*}\|_{0},$

and we conclude that
$\|{\mathbf{y}}\|_{0}=\|\widetilde{{\mathbf{x}}}^{*}\|_{0}$, which together
with (24) complete the argument. ∎

We now show that precise information about the magnitude/strength may make it
possible to (perfectly) recover sources with constant strength

${\mathbf{x}}^{*}=k\sum_{\mathcal{J}}{\mathbf{e}}_{j},\quad k>0,$ (25)

i.e., $x_{j}^{*}=k$ for all $j\in\mathcal{J}=\mathrm{supp}({\mathbf{x}}^{*})$.

######  Corollary 2.1.3.

If there exists a vector ${\mathbf{c}}$ such that (16) and (17) hold, then
${\mathbf{x}}^{*}$ is the unique solution of

$\min_{0\leq{\mathbf{x}}\leq
k}\|{\mathsf{W}}{\mathbf{x}}\|_{1}\quad\textnormal{subject
to}\quad{\mathsf{A}}{\mathbf{x}}={\mathsf{A}}\overbrace{k\sum_{\mathcal{J}}{\mathbf{e}}_{j}}^{={\mathbf{x}}^{*}}.$
(26)

###### Proof.

Problem (26) is the same as problem (18) with a true source ${\mathbf{x}}^{*}$
in constant strength form (25) and with $s=k$. Theorem 2.1 asserts that
${\mathbf{x}}^{*}$ solves (26) and that any other solution ${\mathbf{y}}$ of
this problem must be such that
$\mathrm{supp}({\mathbf{y}})\subseteq\mathrm{supp}({\mathbf{x}}^{*})$. Since
both ${\mathbf{x}}^{*}$ and ${\mathbf{y}}$ solve (26), it follows that
$\|{\mathsf{W}}{\mathbf{y}}\|_{1}=\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}$, and
we find that, see (12),

$\sum_{\mathcal{J}}y_{j}\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}=\sum_{\mathcal{J}}x_{j}^{*}\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}=\sum_{\mathcal{J}}k\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}$

or

$\sum_{\mathcal{J}}(k-y_{j})\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}=0.$

This equality and the constraint $0\leq y_{j}\leq k$, $\forall
j\in\mathcal{J}$, lead to the conclusion that $y_{j}=k$ for all
$j\in\mathcal{J}$ and hence that ${\mathbf{y}}={\mathbf{x}}^{*}$. ∎

In practice one typically does not have exact knowledge about the constant
strength $k$ of the source(s). We now present a result which enables the
recovery, under ideal conditions, of $k$. In what follows, we assume that $k$
is a fixed unknown positive number, and that we want to recover a constant
strength source (25).

First we note that the state equation must have a solution in order for the
basis pursuit problem (8) to be meaningful, and we thus introduce the set

$S=\left\\{s\;|\;\exists{\mathbf{x}}\mbox{ such that
}\|{\mathbf{x}}\|_{\infty}\leq s\mbox{ and
}\quad{\mathsf{A}}{\mathbf{x}}={\mathsf{A}}\,k\sum_{\mathcal{J}}{\mathbf{e}}_{j}\right\\}.$
(27)

######  Corollary 2.1.4.

Assume that there exists a vector ${\mathbf{c}}$ such that (16) and (17) hold
for ${\mathbf{x}}^{*}=k\sum_{\mathcal{J}}{\mathbf{e}}_{j}$,
$\mathcal{J}=\textnormal{supp}({\mathbf{x}}^{*})$. Then the function
$g:S\rightarrow\mathbb{R}$,

$g(s):=\min_{0\leq{\mathbf{x}}\leq
s}\|{\mathsf{W}}{\mathbf{x}}\|_{1}\quad\textnormal{subject
to}\quad{\mathsf{A}}{\mathbf{x}}={\mathsf{A}}\overbrace{k\sum_{\mathcal{J}}{\mathbf{e}}_{j}}^{={\mathbf{x}}^{*}},$

is non-increasing wrt. $s$ and satisfies

$\displaystyle g(s)$
$\displaystyle>\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}\quad\mbox{for
}s\in[0,k)\cap S,$ (28) $\displaystyle g(s)$
$\displaystyle=\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}\quad\mbox{for }s\geq k.$
(29)

###### Proof.

The function $g(s)$ is non-increasing because the size of the constraining set
$\\{{\mathbf{x}}\,|\,0\leq{\mathbf{x}}\leq s\\}$ increases as $s$ increases.

Since $\|{\mathbf{x}}^{*}\|_{\infty}=k$, it follows from Theorem 2.1 that
${\mathbf{x}}^{*}$ solves

$\min_{0\leq{\mathbf{x}}\leq
s}\|{\mathsf{W}}{\mathbf{x}}\|_{1}\quad\textnormal{subject
to}\quad{\mathsf{A}}{\mathbf{x}}={\mathsf{A}}{\mathbf{x}}^{*}$ (30)

when $s\geq k$. Consequently, (29) must hold.

According to Corollary 2.1.3, ${\mathbf{x}}^{*}$ is the unique solution of
(30) when $s=k$. Hence, it follows that (28) must be true. ∎

If the true source(s) has constant strength, we thus obtain the following
scheme for recovering both its size and magnitude $k$:

#### Algorithm

  * •

Compute

$g(s)=\min_{0\leq{\mathbf{x}}\leq
s}\|{\mathsf{W}}{\mathbf{x}}\|_{1}\quad\textnormal{subject
to}\quad{\mathsf{A}}{\mathbf{x}}={\mathbf{b}}$ (31)

many times, varying size of the upper bound $s$.

  * •

The strength $k$ is then the number such that $g(s)=c$ (a constant) for $s\geq
k$ and $g(s)>c$ for $s<k$.

  * •

Corollary 2.1.3 assures the perfect recovery of the true constant strength
source by solving (31) with $s=k$.

We will explore a regularized version of this algorithm in the numerical
experiments section below.

#### Remarks about assumptions (16) and (17)

Our results rely on whether there exists a vector ${\mathbf{c}}$ such that
(16) and (17) hold. This is a subtle question, which we will now discuss in
some detail.

For example, if the cardinality of $\mathcal{J}$ equals two, i.e.,
$\mathcal{J}=\\{j_{1},j_{2}\\}$, then (16) holds with

${\mathbf{c}}=\left(1+\frac{\mathsf{P}{\mathbf{e}}_{j_{1}}}{\|\mathsf{P}{\mathbf{e}}_{j_{1}}\|_{2}}\cdot\frac{\mathsf{P}{\mathbf{e}}_{j_{2}}}{\|\mathsf{P}{\mathbf{e}}_{j_{2}}\|_{2}}\right)^{-1}\left(\frac{\mathsf{P}{\mathbf{e}}_{j_{1}}}{\|\mathsf{P}{\mathbf{e}}_{j_{1}}\|_{2}}+\frac{\mathsf{P}{\mathbf{e}}_{j_{2}}}{\|\mathsf{P}{\mathbf{e}}_{j_{2}}\|_{2}}\right).$

On the other hand, (17) may fail to hold for this particular choice of
${\mathbf{c}}$.

More generally, using the notation
${\mathbf{c}}=\sum_{l=1}^{n}c_{l}{\mathbf{e}}_{l}$, (16) leads to a linear
system with $n$ unknowns $c_{1},\,c_{2},\ldots,\,c_{n}$,

$\sum_{l=1}^{n}{\mathbf{e}}_{l}\cdot\frac{\mathsf{P}{\mathbf{e}}_{j}}{\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}}\,c_{l}-1=0\quad\forall
j\in\mathcal{J},$ (32)

which must have a solution obeying the inequality constraints, see (17),

$\begin{split}&\sum_{l=1}^{n}{\mathbf{e}}_{l}\cdot\frac{\mathsf{P}{\mathbf{e}}_{i}}{\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}}\,c_{l}-1<0\quad\forall
i\in\mathcal{J}^{c}.\end{split}$ (33)

In other words, the solution set of (32) must overlap with the convex domain
specified in (33). This issue can be explored with standard software tools: We
suggest to solve the following linear programming problem (LPP), where
$\delta>0$ is small fixed number,

$\max\sum_{i\in\mathcal{J}^{c}}\epsilon_{i}$ (34)

subject to $\\{c_{l}\\}$ and $\\{\epsilon_{i}\\}$ satisfying

$\begin{split}&\sum_{l=1}^{n}{\mathbf{e}}_{l}\cdot\frac{\mathsf{P}{\mathbf{e}}_{j}}{\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}}\,c_{l}-1=0\quad\forall
j\in\mathcal{J},\\\
&\sum_{l=1}^{n}{\mathbf{e}}_{l}\cdot\frac{\mathsf{P}{\mathbf{e}}_{i}}{\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}}\,c_{l}-1+\delta+\epsilon_{i}=0\quad\forall
i\in\mathcal{J}^{c},\\\ &\epsilon_{i}\geq 0\quad\forall
i\in\mathcal{J}^{c},\end{split}$ (35)

where $\\{\epsilon_{i}\,|\,i\in\mathcal{J}^{c}\\}$ are slack variables. If one
during the solution process of (34)-(35) detects $\\{c_{l}\\}$ and
$\\{\epsilon_{i}\\}$ which satisfy the constraints (35), then one can stop the
LPP program and conclude that the solution set of (32) overlaps with the
convex domain described in (33). That is, there exists a vector ${\mathbf{c}}$
such that (16) and (17) hold. We also note that this approach only yields a
sufficient criteria for asserting that (16) and (17) are satisfied.

Let us also comment that, if the cardinality of the index set $\mathcal{J}$ is
significantly smaller than the number $n$ of columns in the matrix
${\mathsf{A}}\in\mathbb{R}^{m\times n}$, then the solution set of (32) will
typically be relatively large and it becomes more ”likely” that there exists
${\mathbf{c}}$ such that (33) holds. Similarly, if the size of $\mathcal{J}$
increases, then one would expect that the solution set of (32) decreases and
it becomes less ”likely” that (33) is satisfied, i.e., perfect source recovery
with (18) or (26) will ”probably” not be possible. (Here we use ”likely” and
”probably” in an informal manner.)

In the event that there does not exist ${\mathbf{c}}$ such that (32) and (33)
hold, one might wonder whether it is possible to bound the size of the support
of a solution ${\mathbf{y}}$ of (18). That is, to what extent will the size of
the support of the true source
${\mathbf{x}}^{*}=\sum_{\mathcal{J}}x_{j}^{*}{\mathbf{e}}_{j}$ potentially be
increased by solving (18)?

To investigate this issue, assume that $\mathcal{J}^{\prime}$ is an extension
of $\mathcal{J}$ such that properties analogous to (16) and (17) hold: There
exists a vector ${\mathbf{c}}$ such that

$\displaystyle\frac{\mathsf{P}{\mathbf{e}}_{j}}{\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}}\cdot{\mathbf{c}}$
$\displaystyle=1\quad\forall j\in\mathcal{J}^{\prime},$ (36)
$\displaystyle\frac{\mathsf{P}{\mathbf{e}}_{i}}{\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}}\cdot{\mathbf{c}}$
$\displaystyle<1\quad\forall i\in(\mathcal{J}^{\prime})^{c}.$ (37)

Using the same techniques as in the proof of Theorem 2.1, we prove the
following result.

######  Proposition 2.2.

Any solution ${\mathbf{y}}$ of (18) satisfies

$\mathrm{supp}({\mathbf{y}})\subseteq\mathcal{J}^{\prime}$

for any index set
$\mathcal{J}^{\prime}\supseteq\mathcal{J}=\textnormal{supp}({\mathbf{x}}^{*})$
satisfying (36)-(37). Consequently, a solution of (18) must obey

$\mathrm{supp}({\mathbf{y}})\subseteq\overline{\mathcal{J}},$

where

$\overline{\mathcal{J}}=\bigcap\\{\mathcal{J}^{\prime}\;|\;\mathcal{J}^{\prime}\supseteq\mathcal{J}\mbox{
and }\exists\,{\mathbf{c}}\mbox{ such that }\mathcal{J}^{\prime}\mbox{
satisfies \eqref{IC12}-\eqref{IC13}}\\},$

and we assume that $\overline{\mathcal{J}}\neq\emptyset$. This holds for any
$s\in[\|{\mathbf{x}}^{*}\|_{\infty},\infty]$.

###### Proof.

Let ${\mathbf{y}}$ denote a solution of (18). Very similar to the reasoning
leading to (19), we find that

$\sum_{(\mathcal{J}^{\prime})^{c}\cup\mathcal{J}^{\prime}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}\frac{\mathsf{P}{\mathbf{e}}_{i}}{\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}}\cdot{\mathbf{c}}=\sum_{\mathcal{J}}x_{j}^{*}\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}.$

Assuming that there exists $i^{\prime}\in(\mathcal{J}^{\prime})^{c}$ such that
$y_{i}>0$, we can argue as we did in (20)-(21) to conclude that

$\|{\mathsf{W}}{\mathbf{y}}\|_{1}>\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1},$

where ${\mathbf{x}}^{*}=\sum_{\mathcal{J}}x_{j}^{*}{\mathbf{e}}_{j}$. Since
${\mathbf{x}}^{*}$ satisfies the constraints in (18), there can not exist such
$i^{\prime}\in(\mathcal{J}^{\prime})^{c}$. It follows that
$\mathrm{supp}({\mathbf{y}})\subseteq\mathcal{J}^{\prime}$.

∎

As far as we understand, this result is primarily of theoretical interest
because in practice it seems difficult to compute the smallest extension
$\overline{\mathcal{J}}$. We also note that, even though there exists, by
assumption, a vector ${\mathbf{c}}$ for each $\mathcal{J}^{\prime}$ such that
(36) and (37) hold, it does not follow from our argument that
$\overline{\mathcal{J}}$ satisfies these conditions.

###  2.3 Optimality conditions

In the previous subsection we showed that the support of any solution of (18)
is a subset of the support of the true source ${\mathbf{x}}^{*}$ given the
assumptions (16)-(17). We will now show that these assumptions, which we
motivated from a geometrical point view, in fact correspond to the existence
of a KKT/Lagrange multiplier for the optimality system associated with (18).

First, since
${\mathsf{A}}{\mathbf{x}}={\mathsf{A}}{\mathbf{x}}^{*}\Leftrightarrow\mathsf{P}{\mathbf{x}}=\mathsf{P}{\mathbf{x}}^{*}$,
recall that $\mathsf{P}={\mathsf{A}}^{\dagger}{\mathsf{A}}$, we can formulate
the standard KKT-conditions for (18), see, e.g., [24, Theorem 3.34], with
$s=\infty$, as

$\displaystyle{\mathsf{W}}\partial\|{\mathbf{x}}\|_{1}-\mathsf{P}\bm{\lambda}-\bm{\mu}$
$\displaystyle\ni$ $\displaystyle\mathbf{0},$
$\displaystyle\mathsf{P}{\mathbf{x}}-\mathsf{P}{\mathbf{x}}^{*}$
$\displaystyle=$ $\displaystyle\mathbf{0},$
$\displaystyle{\mathbf{x}}\cdot\bm{\mu}$ $\displaystyle=$ $\displaystyle 0,$
$\displaystyle{\mathbf{x}},\bm{\mu}$ $\displaystyle\geq$
$\displaystyle\mathbf{0},$

where $\partial$ denotes the subgradient, and we have used the chain rule and
the fact that the diagonal entries of ${\mathsf{W}}$ are positive:

$\partial_{\mathbf{x}}\|{\mathsf{W}}{\mathbf{x}}\|_{1}={\mathsf{W}}\partial\|{\mathsf{W}}{\mathbf{x}}\|_{1}={\mathsf{W}}\partial\|{\mathbf{x}}\|_{1}.$

Recall that ${\mathbf{x}}^{*}=\sum_{\mathcal{J}}x_{j}^{*}{\mathbf{e}}_{j}$,
${\mathbf{x}}^{*}\geq\mathbf{0}$. Thus, for ${\mathbf{x}}^{*}$ to be a
solution of (18) it follows that there must exist some $\bm{\lambda}^{*}$ and
$\bm{\mu}^{*}$ such that

$[{\mathsf{W}}\partial\|{\mathbf{x}}^{*}\|_{1}]_{i}\ni\begin{cases}[\mathsf{P}\bm{\lambda}^{*}]_{i},\quad&i\in\mathcal{J},\\\
[\mathsf{P}\bm{\lambda}^{*}]_{i}+\mu_{i}^{*},\quad&i\in\mathcal{J}^{c}.\end{cases}$
(38)

(Where $[\mathbf{v}]_{i}$ represents the $i$’th component of the Euclidean
vector $\mathbf{v}$). Here we have used the implication

$\left.\begin{array}[]{l}{\mathbf{x}}^{*}\cdot\bm{\mu}^{*}=0\\\
{\mathbf{x}}^{*},\bm{\mu}^{*}\geq\mathbf{0}\\\
\mathcal{J}=\mathrm{supp}({\mathbf{x}}^{*})\end{array}\right\\}\Longrightarrow\;\;\mu_{i}^{*}=0\;\;\forall
i\in\mathcal{J}.$

If we evaluate the subgradient

$\left[\partial\|{\mathbf{x}}^{*}\|_{1}\right]_{i}=\begin{cases}\\{1\\},\quad&i\in\mathcal{J},\\\
[-1,1],\quad&i\in\mathcal{J}^{c},\end{cases}$

in (38), and combine it with

$[\mathsf{P}\bm{\lambda}^{*}]_{i}=\mathsf{P}\bm{\lambda^{*}}\cdot{\mathbf{e}}_{i}=\bm{\lambda}^{*}\cdot\mathsf{P}{\mathbf{e}}_{i},$

we derive that

$\displaystyle\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}$ $\displaystyle=$
$\displaystyle\mathsf{P}{\mathbf{e}}_{j}\cdot\bm{\lambda}^{*},\quad
j\in\mathcal{J},$ $\displaystyle\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}$
$\displaystyle\geq$
$\displaystyle\mathsf{P}{\mathbf{e}}_{j}\cdot\bm{\lambda}^{*},\quad
j\in\mathcal{J}^{c},$

because
${\mathsf{W}}=\textnormal{diag}\left(\|\mathsf{P}{\mathbf{e}}_{1}\|_{2},\,\|\mathsf{P}{\mathbf{e}}_{2}\|_{2},\,\ldots,\,\|\mathsf{P}{\mathbf{e}}_{n}\|_{2}\right)$,
see (12), and due to the fact that $\bm{\mu}^{*}\geq\mathbf{0}$.

If we have strict complementary slackness, i.e., $\mu^{*}_{i}>0$ whenever
$x^{*}_{i}=0$, we get the stronger result

$\displaystyle\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}$ $\displaystyle=$
$\displaystyle\mathsf{P}{\mathbf{e}}_{j}\cdot\bm{\lambda}^{*},\quad
j\in\mathcal{J},$ $\displaystyle\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}$
$\displaystyle>$
$\displaystyle\mathsf{P}{\mathbf{e}}_{j}\cdot\bm{\lambda}^{*},\quad
j\in\mathcal{J}^{c},$

which are assumptions (16)-(17) with $\mathbf{c}=\bm{\lambda}^{*}$.

###  2.4 Regularized problems

Since $\mathsf{P}={\mathsf{A}}^{\dagger}{\mathsf{A}}$, see (11), the null
spaces of $\mathsf{P}$ and ${\mathsf{A}}$ coincide, and we may write the basis
pursuit problem (8) in the form

$\min_{0\leq{\mathbf{x}}\leq
s}\|{\mathsf{W}}{\mathbf{x}}\|_{1}\quad\textnormal{subject
to}\quad\mathsf{P}{\mathbf{x}}={\mathsf{A}}^{\dagger}{\mathbf{b}}.$ (39)

We also recall that the entries of the regularization matrix ${\mathsf{W}}$
(12) are defined in terms of the projection $\mathsf{P}$. Based on these
observations, and due to mathematical convenience, we will in this subsection
analyse a regularized version of (39) instead of (4):

$\min_{0\leq{\mathbf{x}}\leq
s}\left\\{\frac{1}{2}\|\mathsf{P}{\mathbf{x}}-{\mathsf{A}}^{\dagger}{\mathbf{b}}\|_{2}^{2}+\alpha\sum_{i}w_{i}|x_{i}|\right\\}.$
(40)

More precisely, if
${\mathbf{b}}={\mathbf{b}}^{\dagger}={\mathsf{A}}{\mathbf{x}}^{*}$, to what
extent can we recover ${\mathbf{x}}^{*}$ by solving (40)? Note that, in this
case
${\mathsf{A}}^{\dagger}{\mathbf{b}}^{\dagger}={\mathsf{A}}^{\dagger}{\mathsf{A}}{\mathbf{x}}^{*}=\mathsf{P}{\mathbf{x}}^{*}$.

We remark that we have not managed to prove results similar to those presented
in this subsection for the problem (4). In fact, numerical experiments
indicate that we do not necessarily obtain as good reconstructions with (4) as
with (40) when $\alpha>0$ is not very close to zero.

We now prove a regularized counterpart to Theorem 2.1.

######  Theorem 2.3.

Let ${\mathbf{x}}^{*}=\sum_{\mathcal{J}}x_{j}^{*}{\mathbf{e}}_{j}\geq 0$,
$\mathcal{J}=\mathrm{supp}({\mathbf{x}}^{*})$, and assume that there exist a
vector ${\mathbf{d}}$ and a constant $\gamma<1$ such that

$\displaystyle\frac{\mathsf{P}{\mathbf{e}}_{j}}{\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}}\cdot{\mathbf{d}}$
$\displaystyle=1\quad\forall j\in\mathcal{J},$ (41)
$\displaystyle\frac{\mathsf{P}{\mathbf{e}}_{i}}{\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}}\cdot{\mathbf{d}}$
$\displaystyle<\gamma\quad\forall i\in\mathcal{J}^{c}.$ (42)

Then any

${\mathbf{y}}\in\operatorname*{arg\,min}_{0\leq{\mathbf{x}}\leq
s}\underbrace{\left\\{\frac{1}{2}\left\|\mathsf{P}{\mathbf{x}}-\mathsf{P}\overbrace{\sum_{\mathcal{J}}x_{j}^{*}{\mathbf{e}}_{j}}^{={\mathbf{x}}^{*}}\right\|_{2}^{2}+\alpha\|{\mathsf{W}}{\mathbf{x}}\|_{1}\right\\}}_{:=\mathcal{T}_{\alpha}({\mathbf{x}})}$
(43)

satisfies

$\displaystyle\left|\|{\mathsf{W}}{\mathbf{y}}\|_{1}-\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}\right|$
$\displaystyle\leq c\,\frac{\sqrt{\alpha}}{1-\gamma}$ (44)

and

$\displaystyle\|{\mathsf{W}}{\mathbf{y}}\|_{1,\mathcal{J}^{c}}:=\sum_{\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}$
$\displaystyle\leq c\,\frac{\sqrt{\alpha}}{1-\gamma},$ (45)

provided that $s\in[\|{\mathbf{x}}^{*}\|_{\infty},\infty]$. Here, $c$ is a
constant which is independent of $\alpha$, $\gamma$ and $s$, but depends on
${\mathbf{x}}^{*}$ and ${\mathbf{d}}$.

#### Remark

If there exists a vector ${\mathbf{c}}$ such that (16) and (17) hold, then
there is a vector ${\mathbf{d}}$ and constant $\gamma<1$ such that (41) and
(42) are satisfied because we consider finite dimensional problems:
${\mathbf{d}}={\mathbf{c}}$ and choose $\gamma$ such that

$\displaystyle\max_{i\in\mathcal{J}^{c}}\left\\{\frac{\mathsf{P}{\mathbf{e}}_{i}}{\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}}\cdot{\mathbf{c}}\right\\}<\gamma<1.$

###### Proof of Theorem 2.3.

Assume that ${\mathbf{y}}$ belongs to the $\operatorname*{arg\,min}$ set
defined in (43). Then, because $s\geq\|{\mathbf{x}}^{*}\|_{\infty}$,

$\mathcal{T}_{\alpha}({\mathbf{y}})\leq\mathcal{T}_{\alpha}({\mathbf{x}}^{*})=\alpha\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1},$

which implies that

$\displaystyle\alpha\|{\mathsf{W}}{\mathbf{y}}\|_{1}$
$\displaystyle\leq\alpha\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1},$ (46)
$\displaystyle\|{\bm{\tau}}\|_{2}^{2}$ $\displaystyle\leq
2\alpha\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1},$ (47)

where

${\bm{\tau}}=\mathsf{P}{\mathbf{y}}-\mathsf{P}{\mathbf{x}}^{*}.$

Now,

$\mathsf{P}{\mathbf{y}}\cdot{\mathbf{d}}=\mathsf{P}{\mathbf{x}}^{*}\cdot{\mathbf{d}}+{\bm{\tau}}\cdot{\mathbf{d}},$

which leads to

$\sum_{\mathcal{J}\cup\mathcal{J}^{c}}y_{i}\mathsf{P}{\mathbf{e}}_{i}\cdot{\mathbf{d}}=\sum_{\mathcal{J}}x_{j}^{*}\mathsf{P}{\mathbf{e}}_{j}\cdot{\mathbf{d}}+{\bm{\tau}}\cdot{\mathbf{d}}.$

Combining assumptions (41) and (42) with the constraint $y_{i}\geq 0$ yield
that

$\gamma\sum_{\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}+\sum_{\mathcal{J}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}\geq\sum_{\mathcal{J}}x_{j}^{*}\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}+{\bm{\tau}}\cdot{\mathbf{d}}$

or, also using the fact that
$\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}=\sum_{\mathcal{J}}x_{j}^{*}\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}$
and (46),

$\displaystyle\sum_{\mathcal{J}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}$
$\displaystyle\geq\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}-\gamma\sum_{\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}+{\bm{\tau}}\cdot{\mathbf{d}}$
$\displaystyle\geq\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}-\gamma\left[\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}-\sum_{\mathcal{J}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}\right]+{\bm{\tau}}\cdot{\mathbf{d}}$
$\displaystyle=(1-\gamma)\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}+\gamma\sum_{\mathcal{J}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}+{\bm{\tau}}\cdot{\mathbf{d}},$

which, since $\gamma<1$, implies that

$\sum_{\mathcal{J}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}\geq\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}+\frac{{\bm{\tau}}\cdot{\mathbf{d}}}{1-\gamma}.$
(48)

From (46), keeping in mind that $y_{i}\geq 0$ and invoking (48), it follows
that

$\displaystyle\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}\geq\|{\mathsf{W}}{\mathbf{y}}\|_{1}$
$\displaystyle=\sum_{\mathcal{J}\cup\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}$
$\displaystyle\geq\sum_{\mathcal{J}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}$
$\displaystyle\geq\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}+\frac{{\bm{\tau}}\cdot{\mathbf{d}}}{1-\gamma},$
(49)

and therefore

$\displaystyle\left|\|{\mathsf{W}}{\mathbf{y}}\|_{1}-\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}\right|$
$\displaystyle\leq\left|\frac{{\bm{\tau}}\cdot{\mathbf{d}}}{1-\gamma}\right|$
$\displaystyle\leq\underbrace{\sqrt{2}\|{\mathbf{d}}\|_{2}\sqrt{\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}}}_{:=c}\,\frac{\sqrt{\alpha}}{1-\gamma},$

where we have used the Cauchy-Schwarz inequality and (47). This completes the
proof of (44).

Combing (46) and (48) it follows that

$\displaystyle\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}\geq\|{\mathsf{W}}{\mathbf{y}}\|_{1}$
$\displaystyle=\sum_{\mathcal{J}\cup\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}$
$\displaystyle\geq\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}+\frac{{\bm{\tau}}\cdot{\mathbf{d}}}{1-\gamma}+\sum_{\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2},$

which implies that

$-\frac{{\bm{\tau}}\cdot{\mathbf{d}}}{1-\gamma}\geq\sum_{\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}\geq
0,$

where we have also used the constraint $y_{i}\geq 0$. Hence,

$\displaystyle\sum_{\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}$
$\displaystyle\leq\left|\frac{{\bm{\tau}}\cdot{\mathbf{d}}}{1-\gamma}\right|$
$\displaystyle\leq\sqrt{2}\|{\mathbf{d}}\|_{2}\sqrt{\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}}\,\frac{\sqrt{\alpha}}{1-\gamma},$

which completes the proof of (45). (Note also that (49) implies that
$\frac{{\bm{\tau}}\cdot{\mathbf{d}}}{1-\gamma}\leq 0$.) ∎

The next result contains the regularized version of Corollary 2.1.4. That is,
we consider the possibility for determining the strength $k$ of a source with
constant strength:

${\mathbf{x}}^{*}=k\sum_{\mathcal{J}}{\mathbf{e}}_{j}.$

######  Corollary 2.3.1.

If assumptions (41) and (42) hold with
$\mathcal{J}=\textnormal{supp}({\mathbf{x}}^{*})$, then any

${\mathbf{y}}\in\operatorname*{arg\,min}_{0\leq{\mathbf{x}}\leq
s}\underbrace{\left\\{\frac{1}{2}\left\|\mathsf{P}{\mathbf{x}}-\mathsf{P}\overbrace{k\sum_{\mathcal{J}}{\mathbf{e}}_{j}}^{={\mathbf{x}}^{*}}\right\|_{2}^{2}+\alpha\|{\mathsf{W}}{\mathbf{x}}\|_{1}\right\\}}_{=\mathcal{T}_{\alpha}({\mathbf{x}})}$
(50)

satisfies

$\displaystyle\|{\mathsf{W}}{\mathbf{y}}\|_{1}\geq\left(1+\frac{(1-\gamma)(k-s)}{k\gamma}\right)\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}-C\,\frac{1}{\gamma}\sqrt{s}\sqrt{\alpha}\quad\mbox{for
}s\in[0,k)\cap S,$ (51)
$\displaystyle\left|\|{\mathsf{W}}{\mathbf{y}}\|_{1}-\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}\right|\leq
c\,\frac{\sqrt{\alpha}}{1-\gamma}\quad\mbox{for }s\geq k,$ (52)

where the set $S$ is defined in (27) and $c$ and $C$ are constants which are
independent of $\gamma$, $s$ and $\alpha$.

###### Proof.

Inequality (52) is the same as inequality (44), and we will thus focus on
(51). Let

$s\in[0,k)\cap S$

and assume that ${\mathbf{y}}$ belongs to the $\operatorname*{arg\,min}$ set
in (50). Then, according to the definition (27) of $S$, there exists
${\mathbf{z}}$, with $\|{\mathbf{z}}\|_{\infty}\leq s$, such that
${\mathsf{A}}{\mathbf{z}}={\mathsf{A}}{\mathbf{x}}^{*}$, which implies that
$\mathsf{P}{\mathbf{z}}=\mathsf{P}{\mathbf{x}}^{*}$. It follows that

$\mathcal{T}_{\alpha}({\mathbf{y}})\leq\mathcal{T}_{\alpha}({\mathbf{z}})=\alpha\|{\mathsf{W}}{\mathbf{z}}\|_{1}\leq\alpha
s\|{\mathsf{W}}\mathbf{1}\|_{1},$

which implies that

$\displaystyle\alpha\|{\mathsf{W}}{\mathbf{y}}\|_{1}$ $\displaystyle\leq\alpha
s\|{\mathsf{W}}\mathbf{1}\|_{1},$ (53) $\displaystyle\|{\bm{\tau}}\|_{2}^{2}$
$\displaystyle\leq 2\alpha s\|{\mathsf{W}}\mathbf{1}\|_{1},$ (54)

where

${\bm{\tau}}=\mathsf{P}{\mathbf{y}}-\mathsf{P}{\mathbf{x}}^{*}$

and $\mathbf{1}$ denotes the vector with all components equal to $1$.

Next,

$\mathsf{P}{\mathbf{y}}\cdot{\mathbf{d}}=\mathsf{P}{\mathbf{x}}^{*}\cdot{\mathbf{d}}+{\bm{\tau}}\cdot{\mathbf{d}}$

and it follows that

$\sum_{\mathcal{J}\cup\mathcal{J}^{c}}y_{i}\mathsf{P}{\mathbf{e}}_{i}\cdot{\mathbf{d}}=k\sum_{\mathcal{J}}\mathsf{P}{\mathbf{e}}_{j}\cdot{\mathbf{d}}+{\bm{\tau}}\cdot{\mathbf{d}},$

which we combine with assumptions (41) and (42) and the constraint $y_{i}\geq
0$ to find that

$\gamma\sum_{\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}+\sum_{\mathcal{J}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}\geq
k\sum_{\mathcal{J}}\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}+{\bm{\tau}}\cdot{\mathbf{d}}.$
(55)

It thus follows that

$\sum_{\mathcal{J}\cup\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}\geq(1-\gamma)\sum_{\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}+k\sum_{\mathcal{J}}\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}+{\bm{\tau}}\cdot{\mathbf{d}},$

or

$\|{\mathsf{W}}{\mathbf{y}}\|_{1}\geq(1-\gamma)\sum_{\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}+\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}+{\bm{\tau}}\cdot{\mathbf{d}}$
(56)

because
$\|{\mathsf{W}}{\mathbf{y}}\|_{1}=\sum_{\mathcal{J}\cup\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}$
and
$\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}=k\sum_{\mathcal{J}}\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}$.

Employing the constraint $0\leq{\mathbf{y}}\leq s$ and (55) we find that

$\displaystyle\gamma\sum_{\mathcal{J}^{c}}y_{i}\|\mathsf{P}{\mathbf{e}}_{i}\|_{2}$
$\displaystyle\geq(k-s)\sum_{\mathcal{J}}\|\mathsf{P}{\mathbf{e}}_{j}\|_{2}+{\bm{\tau}}\cdot{\mathbf{d}}$
$\displaystyle=\frac{k-s}{k}\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}+{\bm{\tau}}\cdot{\mathbf{d}},$

which combined with (56) yields

$\displaystyle\|{\mathsf{W}}{\mathbf{y}}\|_{1}$
$\displaystyle\geq(1-\gamma)\frac{(k-s)}{k\gamma}\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}+\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}+\frac{1}{\gamma}{\bm{\tau}}\cdot{\mathbf{d}}$
$\displaystyle\geq\left(1+\frac{(1-\gamma)(k-s)}{k\gamma}\right)\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1}-\frac{1}{\gamma}\sqrt{\alpha}\sqrt{s}\underbrace{\sqrt{2}\sqrt{\|{\mathsf{W}}\mathbf{1}\|_{1}}\|{\mathbf{d}}\|_{2}}_{=C},$

where we have also used (54). This finishes the proof. ∎

##  3 Numerical experiments

We discretized the boundary value problem (2) with the finite element method,
using first-order Lagrange elements, which yields the stiffness matrix
$\mathsf{L}_{\epsilon}$. The forward matrix associated with (1)-(2) then
becomes

${\mathsf{A}}=\mathsf{M}_{\partial}^{1/2}\mathsf{L}_{\epsilon}^{-1}\mathsf{M},$

where $\mathsf{M}$ and $\mathsf{M}_{\partial}$ denote the standard and
”boundary” mass matrices, respectively. With this approach, the local ”hat”
finite element functions $\\{\psi_{i}\\}$ are directly associated with the
standard unit basis vectors $\\{{\mathbf{e}}_{i}\\}$.

Motivated by the discussion presented at the beginning of subsection 2.4, we
performed numerical experiments with

$\min_{0\leq{\mathbf{x}}\leq
s}\left\\{\frac{1}{2}\|\mathsf{P}{\mathbf{x}}-{\mathsf{A}}^{\dagger}{\mathbf{b}}\|_{2}^{2}+\alpha\sum_{i}w_{i}|x_{i}|\right\\}$
(57)

instead of using (4). Here, the matrices
$\mathsf{P}={\mathsf{A}}^{\dagger}{\mathsf{A}}$ and ${\mathsf{W}}$ are defined
in equations (11) and (12), respectively, and the size of the regularization
parameter $\alpha$ is specified in the captions of the figures below. As is
mentioned at the end of subsection 2.1, truncated SVD, using $20$ singular
values, was employed to obtain a computationally reliable approximation of
${\mathsf{A}}^{\dagger}$.

We numerically tested whether one can recover ${\mathbf{x}}^{*}$ by solving
(57) when ${\mathbf{b}}={\mathbf{b}}^{\dagger}={\mathsf{A}}{\mathbf{x}}^{*}$.
If not stated otherwise, the true source always had constant unit strength:

${\mathbf{x}}^{*}=\sum_{\mathcal{J}}{\mathbf{e}}_{j}.$

Since $\mathsf{P}$ has a nontrivial null space and the regularization term in
(57) is not strictly convex, (57) can potentially have several solutions. This
may be ”handled” by adding a small amount of quadratic regularization.
However, we did not do that, and instead we present the outcome of applying
the Alternating Direction Method of Multipliers (ADMM) to (57): Using the zero
initial guess and performing $5000$ iterations. We employed Python, in
combination with the FEniCS software tool, to realize the numerical
experiments.

According to corollaries 2.1.1 and 2.1.3, solving the basis pursuit problem
(18) can, under certain circumstances, lead to exact recovery of the true
source(s), even though ${\mathsf{A}}$ has a large null space. In order to
illuminate these results, we committed an ”inverse crime” in our two first
examples: The same coarse $17\times 17$ grid was used to represent the
source(s) in both the forward and inverse simulations, and the state $u$ was
discretized in terms of a $33\times 33$ mesh. With this setup, the true source
can be exactly represented on the grid used in the inverse computations. This
leads to unrealistically accurate recoveries in the sense that, in addition to
finding the position, size and magnitude, also the exact shape of the
source(s) is recovered – as we will see in examples 1 and 2 below. On the
other had, let us emphasize that, even when inverse crimes are committed, the
recovery of ${\mathbf{x}}^{*}$ is not trivial because ${\mathsf{A}}$ has a
large null space: In example 1 we show that unweighted sparsity regularization
fails to provide adequate results.

In examples 3-5 we avoided inverse crimes: The synthetic boundary data
${\mathbf{b}}$ was generated using a grid with $97\times 97$ nodes for both
the source and the state, whereas a $49\times 49$ grid was used for these
quantities in the inverse solution procedure.

The unit square $\Omega=(0,1)\times(0,1)$ was used as domain, and in all the
cases we ran experiments with both $\epsilon=-1$ and $\epsilon=1$ in (2),
i.e., with both Helmholtz’ equation and the screened Poisson equation. The
results obtained with these two state equations were visually
indistinguishable, and we therefore, in each example, only present results
obtained for one of these model problems: In example 1 and examples 2-5 we
present results for $\epsilon=1$ and $\epsilon=-1$, respectively.

### Example 1: Point sources

In the first example we attempted to locate five single/point sources, see
panel (a) in figure 1. Recall that, if assumptions (16)-(17) hold, then
Theorem 2.1 guarantees that the solution of the basis pursuit problem (18) has
support contained in $\mathcal{J}=\textnormal{supp}({\mathbf{x}}^{*})$ – the
indices associated with the true point sources. Panel (b) in figure 1 shows
the solution of (57) computed with no upper bound, i.e., $s=\infty$. We
observe that all the five point sources, with unit strength, are recovered.
This suggests that ${\mathsf{A}}$ is injective on $\\{{\mathbf{z}}\ |\
\textnormal{supp}({\mathbf{z}})\subset\mathcal{J}\textnormal{ and
}0\leq{\mathbf{z}}\\}$, see Corollary 2.1.1. Indeed, figure 2 shows that also
sources with varying magnitude are perfectly recovered.

(a) True sources.

(b) Inverse solution computed with the constraint $0\leq{\mathbf{x}}$.

(c) Inverse solution computed with the constraint $0\leq{\mathbf{x}}$ and
${\mathsf{W}}=\mathsf{I}$, i.e., unweighted regularization

Figure 1: Example 1. Comparison of the true point sources and the inverse
solutions computed with weighted and unweighted regularization. The
regularization parameter was $\alpha=10^{-4}$.

It was mentioned in the introduction to this section that we committed an
inverse crime in this example. However, we emphasis that the the recovery of
the present point sources is not trivial, which we illustrate by displaying
the solution of

$\min_{0\leq{\mathbf{x}}}\left\\{\frac{1}{2}\|\mathsf{P}{\mathbf{x}}-{\mathsf{A}}^{\dagger}{\mathbf{b}}\|_{2}^{2}+\alpha\sum_{i}|x_{i}|\right\\}$

in panel (c) in figure 1. That is, we ran the same test without weights in the
regularization term. We observe that this procedure fails to recover the true
point sources. Also, since the range of the function plotted in panel (c) is a
subset of $[0,1]$, it would not help to add an upper constraint, i.e, to
employ the box constraint $0\leq{\mathbf{x}}\leq 1$, where we recall that the
strength of the true sources is $1$.

(a) True sources.

(b) Inverse solution computed with the constraint $0\leq{\mathbf{x}}$.

Figure 2: Example 1. Comparison of the true point sources with varying
magnitude and the inverse solution. The regularization parameter was
$\alpha=10^{-4}$. Note that the same grid was used in the forward and inverse
simulations.

### Example 2: Three rectangular sources

Our second example concerns the recovery of three composite sources. More
specifically, we have three separate rectangles, each consisting of several
cells, see panel (a) in figure 3. The right panel of this figure shows that,
when only using the lower constraint ${\mathbf{x}}\geq 0$, our method is
successful in localizing the three sources, although with underestimated
support and overestimated magnitudes. This is in accordance with Theorem 2.1:
The support of the inverse solution is contained in the support of the true
sources. However, if the matrix ${\mathsf{A}}$ is not injective on the set
$\\{{\mathbf{z}}\;|\;\textnormal{supp}({\mathbf{z}})\subset\mathcal{J}\textnormal{
and }0\leq{\mathbf{z}}\leq s\\}$, we can not guarantee exact/perfect recovery.
For $s=\infty$, the matrix ${\mathsf{A}}$ does not appear to be injective on
this set.

To enforce injectivity, we can set the upper box constraint value $s$ equal to
the strength $k$ of the true sources, see Corollary 2.1.3. Nevertheless, in
practise $k$ is typically unknown and it must somehow be estimated. We will
now address this issue.

(a) True sources.

(b) Inverse solution computed with the constraint $0\leq{\mathbf{x}}$.

Figure 3: Example 2. Comparison of the true rectangular sources and the
inverse solution. The regularization parameter was $\alpha=10^{-4}$. Note that
the same grid was used in the forward and inverse simulations.

#### Determining the strength

In connection with Corollary 2.1.4 we presented an algorithm for identifying
the magnitude $k$ of a constant strength source. Analogously, for the
regularized problem (50), one can use inequalities (51)-(52) to approximately
determine $k$:

We can combine (51)-(52) to deduce that a solution ${\mathbf{y}}$ of (57) must
satisfy

$\|{\mathsf{W}}{\mathbf{y}}\|_{1}\sim\begin{cases}\left(1+\frac{(1-\gamma)(k-s)}{k\gamma}\right)\|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1},\
&s\in[0,k)\cap S,\\\ \|{\mathsf{W}}{\mathbf{x}}^{*}\|_{1},\ &s\geq
k,\end{cases}$

provided that $\alpha>0$ is rather small. Thus, we expect that the
$\|{\mathsf{W}}\cdot\|_{1}$-norm of the inverse solution ${\mathbf{y}}$ is
(relatively) constant for $s>k$. Furthermore, for $s<k$, the size of
$\|{\mathsf{W}}\cdot\|_{1}$ should decrease approximately linearly as a
function of $s$. Hence, when plotting $\|{\mathsf{W}}{\mathbf{y}}\|_{1}$ as a
function of the upper bound $s$, we should see an ”L-shaped” curve, which is
indeed what figure 4 displays. Thus, we select the upper box constraint to be
the s-value associated with the vertex of this ”L-curve”.

With this strategy we obtained $s=1$ for the present model problem, which is
indeed the strength of the true sources. The result is displayed in figure
5c). All the three rectangles are exactly identified.

Figure 5 also displays the inverse solutions computed with various other
values of $s$. In panels (a) and (b), $s>1=k$ and, as expected, the support is
a subset of $\mathcal{J}=\textnormal{supp}({\mathbf{x}}^{*})$, see Theorem
2.1. In panels (d)-(f), $s<1=k$, which renders the true source infeasible. For
this particular example, an erroneous choice of $s$ is rather ”forgiving” in
the sense that the support of the inverse solution is either contained inside
the support of the true source (when $s$ is chosen too large) or contained in
a region surrounding the support of the true source (when $s$ is chosen too
small).

Figure 4: Example 2. Plot of the $\|{\mathsf{W}}\cdot\|_{1}$-norm of the
solution ${\mathbf{y}}={\mathbf{y}}(s)$ of (57) for different values of $s$.
Here, $\alpha=10^{-4}.$

### Example 3: The square, the circle and the triangle

In the remaining examples we avoided inverse crimes: A finer grid was employed
in the forward simulations than in the inverse solution procedures. This
yields a more realistic setup, but our analysis is not applicable and exact
recovery is not possible.

Using this setup, can we still recover the size and position of the sources,
but not necessarily the shape? To illuminate the performance of the proposed
method under these circumstances, we considered three sources with different
shape: a square, a rectangle and a circle, see figure 6(a).

The solution computed with only a lower box constraint ${\mathbf{x}}\geq 0$ is
shown in panel (b), and we observe that the support of the inverse solution is
a subset of the support of the true sources. When we also apply the upper box
constraint, then the size of the sources are roughly recovered, but the exact
geometry is ”lost”, as seen in panel (c).

Note that the upper box constraint $s$ equaled $1.2$ in this example. This
value was determined manually by the ”L-curve” procedure described in the
previous example. More precisely, panel (a) in figure 7 shows the
$\|{\mathsf{W}}\cdot\|_{1}$-norm of the inverse solution
${\mathbf{y}}={\mathbf{y}}(s)$ as a function of $s$, and panel (b) displays
the corresponding finite difference approximation of the ”derivative” of this
curve. Compared with the idealized example considered above, cf. figure 4, the
task of selecting an appropriate value for $s$ is now much more challenging.
(Similar observations were made in connection with examples 4 and 5 below.)
The value $1.2$, for $s$, was chosen because the dotted-curve in panel (b) in
figure 7 has a significant ”jump” at this location.

(a) s = 1.4

(b) s = 1.2

(c) s = 1.0

(d) s = 0.8

(e) s = 0.6

(f) s = 0.4

Figure 5: Example 2. Inverse solutions computed with different values for the
upper box constraint $s$. The regularization parameter was $\alpha=10^{-4}$,
and the true rectangular sources are displayed in figure 3(a). The same grid
was used in the forward and inverse simulations.

(a) True sources.

(b) Inverse solution computed with the constraint $0\leq{\mathbf{x}}$.

(c) Inverse solution computed with the box constraint $0\leq{\mathbf{x}}\leq
1.2$.

Figure 6: Example 3. Comparison of the true sources and the inverse solutions.
The regularization parameter was $\alpha=10^{-4}$.

(a) $\|{\mathsf{W}}{\mathbf{y}}(s)\|_{1}$

(b)
$\frac{\|{\mathsf{W}}{\mathbf{y}}(s)\|_{1}-\|{\mathsf{W}}{\mathbf{y}}(s^{\prime})\|_{1}}{s-s^{\prime}}$

Figure 7: Example 3. The $\|{\mathsf{W}}\cdot\|_{1}$-norm of the solution
${\mathbf{y}}={\mathbf{y}}(s)$ of (57), as a function of $s$, and the
associated difference approximation of its ”derivative”. Here,
$\alpha=10^{-4}.$

### Example 4: A horseshoe-shaped source

In examples 4 and 5 we added noise to the data. More specifically, we computed

${\mathbf{b}}={\mathsf{A}}{\mathbf{x}}^{*}+\tau\bm{\rho},$

where $\tau$ is a scalar and $\bm{\rho}$ is a vector containing normally
distributed numbers with zero mean and standard deviation $1$. We refer to the
ratio

$\frac{\tau}{\max\\{{\mathbf{b}}\\}-\min\\{{\mathbf{b}}\\}}$

as the noise level (in percentage). Furthermore, we selected the
regularization parameter $\alpha$ based on Morozov’s discrepancy principle.

The first example with noisy data concerns a horseshoe-shaped source. A
precise recovery of this geometry appears to be more challenging than the
cases considered above because the source now ”partially encloses” a piece of
the domain, see panel (a) in figure 8. This suggests that it might be more
difficult to determine the boundary of the ”inner” part of the horseshoe.

Panels (b)-(d) in figure 8 show that using only the lower constraint
${\mathbf{x}}\geq 0$ results in a highly overestimated magnitude.
Nevertheless, without noise, the support of the solution is a subset of the
support of the true source. When noise is present, the branches (vertical
bars) of the horseshoe is less precisely recovered.

Figure 9 contains results computed with the box constraint
$0\leq{\mathbf{x}}\leq 1$, and we thus used the exact/true upper constraint
$s=1$. In panel (b), the case without noise, we get a rather good
approximation of the horseshoe. For the problem with 1% noise, the branches
are ”bent outwards”, whereas they are almost completely ”gone” for case with
5% noise.

(a) True source.

(b) 0% noise ($\alpha=10^{-4}$).

(c) 1% noise ($\alpha=0.01$).

(d) 5% noise ($\alpha=0.15$).

Figure 8: Example 4. Comparison of the true source and the inverse solutions
computed with the constraint $0\leq{\mathbf{x}}$ for different noise levels.

(a) True source.

(b) 0% noise ($\alpha=10^{-4}$).

(c) 1% noise ($\alpha=0.01$).

(d) 5% noise ($\alpha=0.15$).

Figure 9: Example 4. Comparison of the true source and the inverse solutions
computed with the constraint $0\leq{\mathbf{x}}\leq 1$ for different noise
levels.

### Example 5: A hollow rectangle

In the final example we tried to recover a source shaped like a hollow
rectangle (picture frame), see figure 10(a). This source was deliberately
constructed such that the recovery becomes very challenging. That is, it seems
unrealistic to detect the hole/cavity inside the rectangle when only boundary
data is available.

Indeed, panels (b)-(d) in figure 10 show that our methodology fails to recover
the hole inside the rectangle. On the other hand, the algorithm is quite
successful in determining the ”outer” boundary of the true source for the
cases without noise and with 1% noise. When the data contains 5% noise, both
the position and the size are reasonably approximated, but the shape is too
”circular”.

(a) True source.

(b) 0% noise ($\alpha=10^{-4}$).

(c) 1% noise ($\alpha=0.05$).

(d) 5% noise ($\alpha=0.2$).

Figure 10: Example 5. Comparison of the true source and the inverse solutions
computed with the box constraint $0\leq{\mathbf{x}}\leq 1$ for different noise
levels.

##  4 Concluding remarks

The main purpose of this paper has been to address the non-uniqueness issue
arising in connection with some inverse problems. We have explored this matter
for the task of identifying the source term in an elliptic PDE from boundary
data. Despite the presence of a large null space, we have shown that the use
of weighted sparsity regularization and box constraints can lead to rather
accurate recovery. In fact, in an idealized setting, exact reconstruction
might be possible. We demonstrated this behavior for two test problems, where
we, for the sake of exemplifying the exact recovery, had to commit the
”inverse crime” of using identical grids for the forward and inverse
computations. Under less ideal conditions, i.e., without ”inverse crimes”, we
still observe that the position and roughly the size of the sources can be
computed. That is, the shape can not be perfectly reconstructed – which is as
one would expect.

For the sake of generality, we developed the methodology and the analysis in
terms of Euclidean spaces. Our approach can thus be applied to a broad range
of inverse problems, assuming that one seeks sparse solutions. We presented
mathematical results for both the basis pursuit formulation of the problem and
for its regularized version. It remains to develop an infinite dimensional
counterpart to the theory presented in this paper and to test the performance
of the methodology on real-world data.

This investigation was motivated by the fact that standard (unweighted)
sparsity regularization do not provide adequate results for the source
identification task considered in this paper.

## References

  * [1] B. Abdelaziz, A. El Badia, and A. El Hajj.  Direct algorithms for solving some inverse source problems in 2D elliptic equations.  Inverse Problems, 31(10):105002, 2015.
  * [2] A. Ben Abda, F. Ben Hassen, J. Leblond, and M. Mahjoub.  Sources recovery from boundary data: A model related to electroencephalography.  Mathematical and Computer Modelling, 49:2213–2223, 2009.
  * [3] X. Cheng, R. Gong, and W. Han.  A new Kohn-Vogelius type formulation for inverse source problems.  Inverse Problems and Imaging, 9(4):1051–1067, 2015.
  * [4] A. El Badia and T. Ha-Duong.  An inverse source problem in potential analysis.  Inverse Problems, 16:651–663, 2000.
  * [5] O. L. Elvetun and B. F. Nielsen.  A regularization operator for source identification for elliptic PDEs.  Inverse Problems and Imaging, 15(4):599–618, 2021.
  * [6] O. L. Elvetun and B. F. Nielsen.  Modified Tikhonov regularization for identifying several sources.  International Journal of Numerical Analysis and Modeling, 18(6):740–757, 2021.
  * [7] H. W. Engl, M. Hanke, and A. Neubauer.  Regularization of Inverse Problems.  Kluwer Academic Publishers, 1996.
  * [8] R. Griesmaier, M. Hanke, and T. Raasch.  Inverse source problems for the Helmholtz equation and the windowed Fourier transform.  SIAM Journal on Scientific Computing, 34(3):A1544–A1562, 2012.
  * [9] R. Griesmaier, M. Hanke, and T. Raasch.  Inverse source problems for the Helmholtz equation and the windowed Fourier transform II.  SIAM Journal on Scientific Computing, 35(5):A2188–A2206, 2013\.
  * [10] M. Hanke, L. Harhanen, N. Hyvönen, and E. Schweickert.  Convex source support in three dimensions.  BIT Numerical Mathematics, 52:45–63, 2012.
  * [11] M. Hanke, N. Hyvönen, and S. Reusswig.  Convex source support and its application to electric impedance tomography.  SIAM Journal on Imaging Sciences, 1(4):364–378, 2008.
  * [12] M. Hanke and W. Rundell.  On rational approximation methods for inverse source problems.  Inverse Problems and Imaging, 5(1):185–202, 2011.
  * [13] F. Hettlich and W. Rundell.  Iterative methods for the reconstruction of an inverse potential problem.  Inverse Problems, 12:251–266, 1996.
  * [14] M. Hinze, B. Hofmann, and T. N. T. Quyen.  A regularization approach for an inverse source problem in elliptic systems from single Cauchy data.  Numerical Functional Analysis and Optimization, 40(9):1080–1112, 2019.
  * [15] M. Ikehata.  Reconstruction of a source domain from the Cauchy data.  Inverse Problems, 15(2):637–645, 1999.
  * [16] M. Ikehata.  Reconstruction of a source domain from the Cauchy data: II. three-dimensional case.  Inverse Problems, 37(12):125004, 2021.
  * [17] V. Isakov.  Inverse Problems for Partial Differential Equations.  Springer-Verlag, 2005.
  * [18] R. Kress and W. Rundell.  Reconstruction of extended sources for the Helmholtz equation.  Inverse Problems, 29(3):035005, 2013.
  * [19] K. Kunisch and X. Pan.  Estimation of interfaces from boundary measurements.  SIAM Journal on Control and Optimization, 32(6):1643–1674, 1994\.
  * [20] S. Kusiak and J. Sylvester.  The scattering support.  Communications on Pure and Applied Mathematics, 56(11):1525–1548, 2003.
  * [21] S. Kusiak and J. Sylvester.  The convex scattering support in a background medium.  SIAM Journal on Mathematical Analysis, 36(4):1142–1158, 2005.
  * [22] R. Potthast, J. Sylvester, and S. Kusiak.  A ’range test’ for determining scatterers with unknown physical properties.  Inverse Problems, 19(3):533–547, 2003.
  * [23] W. Ring.  Identification of a core from boundary data.  SIAM Journal on Applied Mathematics, 55(3):677–706, 1995.
  * [24] A. Ruszczyński.  Nonlinear Optimization.  Princeton University Press, 2006.
  * [25] S. J. Song and J. G. Huang.  Solving an inverse problem from bioluminescence tomography by minimizing an energy-like functional.  Journal of Computational Analysis and Applications, 14:544–558, 2012.
  * [26] X. Wang, Y. Guo, D. Zhang, and H. Liu.  Fourier method for recovering acoustic sources from multi-frequency far-field data.  Inverse Problems, 33(3), 2017.
  * [27] D. Zhang, Y. Guo, J. Li, and H. Liu.  Retrieval of acoustic sources from multi-frequency phaseless data.  Inverse Problems, 34(9), 2018.
  * [28] D. Zhang, Y. Guo, J. Li, and H. Liu.  Locating multiple multipolar acoustic sources using the direct sampling method.  Communications in Computational Physics, 25(5):1328–1356, 2019\.