# A lower bound on the tree-width of graphs with irrelevant vertices111This
version has been created by modifying a preprint submitted to Elsevier;
modifications are minor and were done purely to work around lack of support
for xelatex on arXiv

Isolde Adler<EMAIL_ADDRESS>Philipp Klaus Krause
<EMAIL_ADDRESS>School of Computing, University of Leeds
Albert-Ludwigs-Universität Freiburg

###### Abstract

For their famous algorithm for the disjoint paths problem, Robertson and
Seymour proved that there is a function $f$ such that if the tree-width of a
graph $G$ with $k$ pairs of terminals is at least $f(k)$, then $G$ contains a
solution-irrelevant vertex (Graph Minors. XXII., JCTB 2012). We give a single-
exponential lower bound on $f$. This bound even holds for planar graphs.

###### keywords:

disjoint paths problem , irrelevant vertex , vital linkage , unique linkage ,
planar graph , tree-width

##  1 Introduction

The Disjoint Paths Problem is one of the famous classical problems in the area
of graph algorithms. Given a graph $G$, and $k$ pairs of terminals,
$(s_{1},t_{1}),\ldots,(s_{k},t_{k})$, it asks whether $G$ contains $k$ vertex-
disjoint paths $P_{1},\ldots,P_{k}$ such that $P_{i}$ connects $s_{i}$ to
$t_{i}$, (for $i=1,\ldots,k$). Karp proved that the problem is NP-hard in
general [4] and Lynch proved that it remains NP-hard on planar graphs [6].
Robertson and Seymour showed that it can be solved in time
$g(k)\cdot|V(G)|^{3}$ for some computable function $g$, i. e. the problem is
fixed-parameter tractable (and, in particular, solvable in polynomial time for
fixed $k$). For a recursive step in their algorithm ((10.5) in [11]), they
prove [13] that there is a function $f\colon\mathbb{N}\to\mathbb{N}$ such that
if a graph $G$ with $k$ pairs of terminals has tree-width at least $f(k)$,
then $G$ contains a vertex that is _irrelevant_ to the solution, i. e. $G$
contains a non-terminal vertex $v$ such that $G$ has a solution if and only if
the graph $G-v$ (with the same terminals) has a solution.

In this paper we give a lower bound on $f$, showing that $f(k)\geq 2^{k}$,
even for planar graphs. For this we construct a family of planar input graphs
$(G_{k})_{k\geq 2}$, each with $k$ pairs of terminals, such that the tree-
width of $G_{k}$ is $2^{k}-1$, and every member of the family has a unique
solution to the Disjoint Paths Problem, where the paths of the solution use
all vertices of the graph. Hence no vertex of $G_{k}$ is irrelevant. As a
corollary, we obtain a lower bound of $2^{k}-1$ on the tree-width of graphs
having _vital linkages_ (also called _unique linkages_) [12] with $k$
components.222This result appeared in the last section of a conference paper
[1]. While the main focus of the paper [1] was a single exponential upper
bound on $f$ on planar graphs, it only sketches the lower bound. Here we
provide the full proof of the lower bound. A longer proof of the lower bound
can also be found in the thesis [5]. Our result contrasts the polynomial upper
bound in a related topological setting [7], where two systems of curves are
untangled on a sphere with holes.

For planar graphs, an upper bound of $f(k)\leq 72\sqrt{2}k^{\frac{3}{2}}\cdot
2^{k}$ was given in [1]. An elementary proof for a bound of $f(k)\leq(72k\cdot
2^{k}-72\cdot 2^{k}+18)\lceil\sqrt{2k+1}\rceil$ was provided later [5] as well
as a slightly improved bound of $f(k)\leq 26k\cdot 2^{\frac{3}{2}}\cdot 2^{k}$
requiring a slightly more involved proof [2]. Our lower bound shows that this
is asymptotically optimal. Recently, an explicit upper bound on $f$ on graphs
of bounded genus [3] was found, then refined into one that is single
exponential in $k$ and the genus [8]. The exact order of growth of $f$ on
general graphs is still unknown.

##  2 Preliminaries

Let $\mathbb{N}$ denote the set of all non-negative integers. For
$k\in\mathbb{N}$, we let $[k]\mathrel{\mathop{:}}=\\{1,\ldots,k\\}$. For a set
$S$ we let $2^{S}$ denote the power set of $S$. A _graph_ $G=(V,E)$ is a pair
of a set of _vertices_ $V$ and a set of _edges_ $E\subseteq\\{e\ |\ e\in
2^{V},|e|=2\\}$, i. e. graphs are undirected and simple. For an edge
$e=\\{x,y\\}$, the vertices $x$ and $y$ are called _endpoints_ of the edge
$e$, and the edge is said to be between its endpoints. For a graph $G=(V,E)$
let $V(G)\mathrel{\mathop{:}}=V$ and $E(G)\mathrel{\mathop{:}}=E$. Let $H$ and
$G$ be graphs. The graph $H$ is a _subgraph_ of $G$ (denoted by $H\subseteq
G$), if $V(H)\subseteq V(G)$ and $E(H)\subseteq E(G)$. For a set $X\subseteq
V(G)$, the subgraph of $G$ _induced by_ $X$ is the graph
$G[X]\mathrel{\mathop{:}}=(X,\\{e\in E(G)\ |\ e\subseteq X\\})$ and we let
$G-v\mathrel{\mathop{:}}=G[V(G)\smallsetminus\\{v\\}]$.

A _path_ $P$ in a graph $G=(V,E)$ is a sequence $n_{0},\ldots,n_{k}\in V$ of
pairwise distinct vertices of $G$, such that for every
$i\in\\{0,\ldots,k-1\\}$ there is an edge $\\{n_{i},n_{i+1}\\}\in E$. The
vertices $n_{0}$ and $n_{k}$ are called _endpoints_ of $P$. The path $P$ is
called a path _from_ $n_{0}$ _to_ $n_{k}$ (i. e. paths are _simple_). We
sometimes identify the path $P$ in $G$ with the subgraph
$(\\{n_{0},\ldots,n_{k}\\},\\{\\{n_{0},n_{1}\\},\ldots,\\{n_{k-1},n_{k}\\}\\})$
of $G$. A graph $G$ is called _connected_ , if it has at least one vertex and
for any two vertices $x,y\in V(G)$, there is a path from $x$ to $y$ in $G$.
The inclusion-maximal connected subgraphs of a graph are called _connected
components_ of the graph. For $A,B\subseteq V(G)$, a set $S\subseteq V(G)$
separates $A$ from $B$, if there is no path from a vertex in $A$ to a vertex
in $B$ in the subgraph of $G$ induced by $V(G)\smallsetminus S$. A _tree_ is a
non-empty graph $T$, such that for any two vertices $x,y\in V(T)$ there is
exactly one path from $x$ to $y$ in $T$.

Figure 1: $(3\times 3)$-grid

Let $m,n\in\mathbb{N}\smallsetminus\\{0\\}$. The $(m\times n$)-_grid_ is a
graph $H=(V,E)$ with $V\mathrel{\mathop{:}}=[m]\times[n]$ and
$E\mathrel{\mathop{:}}=\\{\\{(y,x),(w,z)\\}\ |\ (y,x)\in V,(w,z)\in
V,|x-z|+|y-w|=1\\}$. In case of a square grid where $m=n$, we say that $n$ is
the _size_ of the grid. An edge $\\{(y,x),(w,z)\\}$ in the grid is called
_horizontal_ , if $y=w$, and _vertical_ , if $x=z$. See Figure 1 for the
$(3\times 3)$-grid.

A _drawing_ of a graph $G$ is a representation of $G$ in the Euclidean plane
$\mathbb{R}^{2}$, where vertices are represented by distinct points of
$\mathbb{R}^{2}$ and edges by simple curves joining the points that correspond
to their endpoints, such that the interior of every curve representing an edge
does not contain points representing vertices. A _planar drawing_ (or
_embedding_) is a drawing, where the interiors of any two curves representing
distinct edges of $G$ are disjoint. A graph $G$ is _planar_ , if $G$ has a
planar drawing (See [10] for more details on planar graphs). A _plane graph_
is a planar graph $G$ together with a fixed embedding of $G$ in
$\mathbb{R}^{2}$. We will identify a plane graph with its image in
$\mathbb{R}^{2}$. Once we have fixed the embedding, we will also identify a
planar graph with its image in $\mathbb{R}^{2}$.

######  Definition 1 (Disjoint Paths Problem ($\operatorname{DPP}$)).

Given a graph $G$ and $k$ pairs of _terminals_ $(s_{1},t_{1})\in
V(G)^{2},\ldots,(s_{k},t_{k})\in V(G)^{2}$, the Disjoint Paths Problem is the
problem of deciding whether $G$ contains $k$ vertex-disjoint paths
$P_{1},\ldots,P_{k}$ such that $P_{i}$ connects $s_{i}$ to $t_{i}$ (for
$i\in[k]$). If such paths $P_{1},\ldots,P_{k}$ exist, we refer to them as a
_solution_. We denote an _instance_ of $\operatorname{DPP}$ by
$G,(s_{1},t_{1}),(s_{2},t_{2}),\ldots,(s_{k},t_{k})$.

Let $G,(s_{1},t_{1}),\ldots,(s_{k},t_{k})$ be an instance of
$\operatorname{DPP}$. A non-terminal vertex $v\in V(G)$ is _irrelevant_ , if
$G,(s_{1},t_{1}),\ldots,(s_{k},t_{k})$ has a solution if and only if
$G-v,(s_{1},t_{1}),\ldots,(s_{k},t_{k})$ has a solution.

A _tree-decomposition_ of a graph $G$ is a pair $(T,\chi)$, consisting of a
tree $T$ and a mapping $\chi\colon V(T)\to 2^{V(G)}$, such that for each $v\in
V(G)$ there exists $t\in V(T)$ with $v\in\chi(t)$, for each edge $e\in E(G)$
there exists a vertex $t\in V(T)$ with $e\subseteq\chi(t)$, and for each $v\in
V(G)$ the set $\\{t\in V(T)\mid v\in\chi(t)\\}$ is connected in $T$. The
_width_ of a tree-decomposition $(T,\chi)$ is

$\operatorname{w}(T,\chi)\mathrel{\mathop{:}}=\max\left\\{\left|\chi(t)\right|-1\
\middle|\ t\in V(T)\right\\}.$

If $T$ is a path, $(T,\chi)$ is also called a _path-decomposition_. The _tree-
width of $G$_ is

$\operatorname{tw}(G)\mathrel{\mathop{:}}=\min\left\\{\operatorname{w}(T,\chi)\
\middle|\ (T,\chi)\text{ is a tree-decomposition of }G\right\\}.$

The _path-width of $G$_ is

$\operatorname{pw}(G)\mathrel{\mathop{:}}=\min\left\\{\operatorname{w}(T,\chi)\
\middle|\ (T,\chi)\text{ is a path-decomposition of }G\right\\}.$

Obviously, every graph $G$ satisfies
$\operatorname{pw}(G)\geq\operatorname{tw}(G)$. Every tree has tree-width at
most $1$ and every path has path-width at most $1$. It is well known that the
$(n\times n$)-grid has both tree-width and path-width $n$. Moreover, if
$H\subseteq G$, then $\operatorname{tw}(H)\leq\operatorname{tw}(G)$ and
$\operatorname{pw}(H)\leq\operatorname{pw}(G)$.

######  Theorem 1 (Robertson and Seymour [13]).

There is a function $f\colon\mathbb{N}\to\mathbb{N}$ such that if
$\operatorname{tw}(G)\geq f(k)$, then $G,(s_{1},t_{1}),\ldots,(s_{k},t_{k})$
has an irrelevant vertex (for any choice of terminals
$(s_{1},t_{1}),\ldots,(s_{k},t_{k})$ in $G$).

A _linkage_ in a graph $G$ is a subgraph $L\subseteq G$, such that each
connected component of $L$ is a path. The _endpoints_ of a linkage $L$ are the
endpoints of these paths, and the _pattern_ of $L$ is the matching on the
endpoints induced by the paths, i. e. the pattern is the set

$\left\\{\\{s,t\\}\ \middle|\ L\text{ has a connected component that is a path
from $s$ to $t$}\right\\}.$

A linkage $L$ in a graph $G$ is a _vital linkage_ in $G$, if $V(L)=V(G)$ and
there is no other linkage $L^{\prime}\neq L$ in $G$ with the same pattern as
$L$.

######  Theorem 2 (Robertson and Seymour [13]).

There are functions $g,h\colon\mathbb{N}\to\mathbb{N}$ such that if a graph
$G$ has a vital linkage with $k$ components then $\operatorname{tw}(G)\leq
g(k)$ and $\operatorname{pw}(G)\leq h(k)$.

##  3 The lower bound

Our main result is the following.

######  Theorem 3.

Let $f,g,h\colon\mathbb{N}\to\mathbb{N}$ be as in Theorems 1 and 2. Then
$f(k)\geq 2^{k}$, $g(k)\geq 2^{k}-1$, and $h(k)\geq 2^{k}-1$. Moreover, this
holds even if we consider planar graphs only.

Figure 2: $G_{4},(s_{1},t_{1}),(s_{2},t_{2}),(s_{3},t_{3}),(s_{4},t_{4})$ with
solution.

In our proof we construct a family of graphs $G_{k},k\geq 1$, of tree-width
and path-width $\geq 2^{k}-1$, and with a vital linkage with $k$ components.
Figure 2 shows the graph $G_{4}$.

Figure 3: The construction of $G_{4}=G_{4,15}$ from $G_{3,30}$.

######  Definition 2 (The graph $G_{k}$).

Let $k,p\in\mathbb{N}\smallsetminus\\{0\\}$. We inductively define an instance
$G_{k},(s_{1},t_{1}),\ldots,(s_{k},t_{k})$ of $\operatorname{DPP}$ as follows.

The Graph $G_{1,p}$ is the path $x_{1},x_{2},\ldots,x_{p}$ with $p$ vertices,
$s_{1}(G_{1,p})\mathrel{\mathop{:}}=x_{1}$,
$t_{1}(G_{1,p})\mathrel{\mathop{:}}=x_{p}$. The _bottom row_ and the _top row_
of $G_{1,p}$ are the graph $G_{1,p}$ itself.

We define the graph $G_{k+1,p}$ by adding a path $y_{1},y_{2},\ldots,y_{p}$
with $p$ vertices to $G_{k,2p}$ as follows. Let $x_{1},x_{2},\ldots,x_{2p}$ be
the bottom row of $G_{k,2p}$ and let $z_{1},z_{2},\ldots,z_{2p}$ be the top
row of $G_{k,2p}$. Let

$\displaystyle V(G_{k+1,p})\mathrel{\mathop{:}}=V(G_{k,2p})$
$\displaystyle\cup\\{y_{1},y_{2},\ldots,y_{p}\\},$ $\displaystyle
E(G_{k+1,p})\mathrel{\mathop{:}}=E(G_{k,2p})$
$\displaystyle\cup\left\\{\\{y_{i},y_{i+1}\\}\ \middle|\ 1\leq
i<p\right\\}\cup$
$\displaystyle\left\\{\\{y_{i},x_{i}\\},\\{y_{i},x_{2p-i+1}\\}\ \middle|\
1\leq i\leq p\right\\}.$

We set $s_{k+1}(G_{k+1,p})\mathrel{\mathop{:}}=y_{1}$,
$t_{k+1}(G_{k+1,p})\mathrel{\mathop{:}}=y_{p}$ and
$s_{i}(G_{k+1,p})\mathrel{\mathop{:}}=s_{i}(G_{k,p})$,
$t_{i}(G_{k+1,p})\mathrel{\mathop{:}}=t_{i}(G_{k,p})$ for $1\leq i\leq k$. The
_top row_ of $G_{k+1,p}$ is $z_{1},\ldots,z_{p}$ and the _bottom row_ of
$G_{k+1,p}$ is $z_{2p},\ldots,z_{p+1}$.

Let $G_{k}\mathrel{\mathop{:}}=G_{k,2^{k}-1}$. We define the
$\operatorname{DPP}$ instance $G_{k},(s_{1},t_{1}),\ldots,(s_{k},t_{k})$ as
$G_{k},(s_{1}(G_{k}),t_{1}(G_{k})),\ldots,(s_{k}(G_{k}),t_{k}(G_{k}))$.

Figure 3 shows the construction of $G_{4}=G_{4,15}$ from $G_{3,30}$.

######  Remark 1.

By construction, the graph $G_{k}$ contains a
$((2^{k}-1)\times(2^{k}-1))$-grid as a subgraph. The tree-width and path-width
of $G_{k}$ are thus at least $2^{k}-1$.

######  Remark 2.

By construction, the graph $G_{k}$ contains a linkage (because in each step we
add a path linking a new terminal pair).

We will now show that this linkage is vital by considering a topological
version.

######  Definition 3 (Topological DPP).

Given a subset $X$ of the plane and $k$ pairs of terminals $(s_{1},t_{1})\in
X^{2},\ldots,(s_{k},t_{k})\in X^{2}$ the topological Disjoint Paths Problem is
the problem of deciding whether there are $k$ pairwise disjoint curves in $X$,
such that each curve $P_{i}$ is homeomorphic to $[0,1]$ and its ends are
$s_{i}$ and $t_{i}$. If such curves $P_{1},\ldots,P_{k}$ exist, we refer to
them as a _solution_. We denote an _instance_ of the topological Disjoint
Paths Problem by $X,(s_{1},t_{1}),(s_{2},t_{2}),\ldots,(s_{k},t_{k})$.

A _disc-with-edges_ is a subset $X$ of the plane containing a closed disc $D$
such that the connected components of $X\smallsetminus D$, called _edges_ ,
are homeomorphic to open intervals $(0,1)$. We now define a family
$(X_{k})_{k\in\mathbb{N}\smallsetminus\\{0\\}}$ of discs-with-edges together
with terminals. These will be used as instances of the topological
$\operatorname{DPP}$. Figure 4 illustrates the construction.

Figure 4: The construction of $X_{4},(s_{1},t_{1}),\ldots,(s_{4},t_{4})$, for
the topological $\operatorname{DPP}$. Note that $s_{5}$ and $s_{6}$ are only
used to place $E_{4}$ correctly.

######  Definition 4 ($X_{k}$).

Let $D$ be a closed disc in the plane and
$k\in\mathbb{N}\smallsetminus\\{0\\}$. We start by inductively defining points
$s_{k}$, $t_{k}$ on the boundary $\partial D$ of $D$. (These will be used as
terminals and to confine the way the edges are added to $D$.) Let
$s_{1},t_{1}$ be two distinct points on $\partial D$, and let
$C_{1}\mathrel{\mathop{:}}=\partial D\smallsetminus\\{s_{1},t_{1}\\}$. Hence
$C_{1}$ is the union of two curves, each homeomorphic to the open interval
$(0,1)$. Call one of the curves $S_{1}$ and the other $T_{1}$. Assume that
$s_{k}$, $t_{k}$, $C_{k}$, $S_{k}$, and $T_{k}$ are already defined, and
assume that $T_{k}$ is a curve adjacent to $t_{k}$ and $s_{1}$. Place a new
point $s_{k+1}$ on $S_{k}$ and a new point $t_{k+1}$ on $T_{k}$, let
$C_{k+1}\mathrel{\mathop{:}}=C_{k}\smallsetminus\\{s_{k+1},t_{k+1}\\}$, let
$T_{k+1}$ be the component of $C_{k+1}$ adjacent to $t_{k+1}$ and $s_{1}$, and
let $S_{k+1}$ be the component of $C_{k+1}$ adjacent to $t_{k+1}$ and $t_{k}$.

Now let $X_{1}\mathrel{\mathop{:}}=D$ and
$E_{1}\mathrel{\mathop{:}}=\varnothing$. Assume the space $X_{k}$ and the set
$E_{k}$ are already defined. We define $X_{k+1}$ by adding a planar matching
of $2^{k}-1$ edges to $X_{k}$. We call the set of these edges $E_{k+1}$. The
edges are pairwise disjoint and disjoint from $X_{k}$. They are added such
that each end is adjacent to a point on $\partial D$ and no two edges are
adjacent to the same point on $\partial D$. Each edge has one end adjacent to
a point on the component of $C_{k+2}$ between $t_{k}$ and $s_{k+1}$, and the
other end adjacent to a point on the component of $C_{k+2}$ between $t_{k}$
and $s_{k+2}$. Finally, let $X_{k+1}\mathrel{\mathop{:}}=X_{k}\cup E_{k+1}$.

In this way we obtain a family $X_{k}$,
$(s_{1},t_{1}),(s_{2},t_{2}),\ldots,(s_{k},t_{k})$ of instances to the
topological DPP.

Figure 5: A solution of the topological $\operatorname{DPP}$ from Figure 4.

######  Remark 3.

The embedding of $G_{k}$ (as shown in Figure 2 for $G_{4}$) corresponds to the
space $X_{k}$. Thus by Remark 2 the topological $\operatorname{DPP}$ on
$X_{k}$, $(s_{1},t_{1}),(s_{2},t_{2}),\ldots,$ $(s_{k},t_{k})$ has a solution.

For an instance of the topological $\operatorname{DPP}$ on $X_{4}$, this
solution can be seen in Figure 5.

######  Lemma 1.

For $k\in\mathbb{N}\smallsetminus\\{0\\}$ the topological $\operatorname{DPP}$
instance $X_{k},(s_{1},t_{1}),\ldots,(s_{k},t_{k})$ has a unique solution
$P_{1},\ldots,P_{k}$ (up to homeomorphism). The solution uses all edges
$\bigcup_{1\leq i\leq k}E_{i}$.

Figure 6: The solution on $X^{\prime}_{3},(s_{1},t_{1}),\ldots,(s_{3},t_{3})$
induced by the solution on $X_{4},(s_{1},t_{1}),\ldots,(s_{4},t_{4})$.

###### Proof.

For $k=1$ this is true because $E_{1}=\varnothing$. Inductively assume that
the lemma holds for $k$. Let $P_{1},\ldots,P_{k+1}$ be any solution to
$X_{k+1},(s_{1},t_{1}),\ldots$, $(s_{k+1},t_{k+1})$. This solution induces a
solution of the topological $\operatorname{DPP}$ $X_{k},(s_{1},t_{1})$,
$\ldots$, $(s_{k},t_{k})$ as follows. Every edge $e\in E_{k+1}$ together with
the segment of $\partial D$ that connects the ends of $e$ and contains
$t_{k+1}$ bounds a disc $D_{e}$. The space
$X^{\prime}_{k}\mathrel{\mathop{:}}=X_{k+1}\cup\bigcup_{e\in E_{k+1}}D_{e}$ is
homeomorphic to $X_{k}$ and the paths $P_{1},\ldots,P_{k}$ form a solution of
$X^{\prime}_{k},$ $(s_{1},t_{1})$, $\ldots$, $(s_{k},t_{k})$. Figure 6
illustrates this for $k=3$. By induction, this solution is unique up to
homeomorphism and the paths $P_{1},\ldots,P_{k}$ use all edges in
$\bigcup_{1\leq i\leq k}E_{i}$. Let $Q_{1},\ldots,Q_{k}$ be the solution
obtained by embedding the graph $G_{k}$ (cf. Remark 3). By uniqueness, for
each $i\in[k]$, the edges of $\bigcup_{1\leq i\leq k}E_{i}$ used by $P_{i}$
are the same as for $Q_{i}$, and the order of their appearance on $P_{i}$ when
walking from $s_{i}$ to $t_{i}$ is also the same as on $Q_{i}$. Hence the
solution $P_{1},\ldots,P_{k}$ on $X^{\prime}_{k}$ restricted to the closed
disc $D$ of $X^{\prime}_{k}$ is a planar matching of curves (the curves in
$\bigcup_{1\leq i\leq k}P_{i}\smallsetminus\bigcup_{1\leq i\leq k}E_{i}$)
between pairs of points on $\partial D$ (and the same pairs of points are
obtained by restricting $Q_{1},\ldots,Q_{k}$ to $D$). These pairs of points
also have to be matched in $X_{k+1}$.

We now claim that in the solution $P_{1},\ldots,P_{k+1}$ on $X_{k+1}$, each
curve in $\bigcup_{1\leq i\leq k}P_{i}\smallsetminus\bigcup_{1\leq i\leq
k}E_{i}$ uses an edge of $E_{k+1}$. If not, then there is a curve

$p\in\bigcup_{1\leq i\leq k}P_{i}\smallsetminus\bigcup_{1\leq i\leq k}E_{i}$

that avoids all edges in $E_{k+1}$. Since the edges of $\bigcup_{1\leq i\leq
k}E_{i}$ are already used, $p$ is routed within $D$. By construction of
$X_{k+1}$ and the fact that all edges of $\bigcup_{1\leq i\leq k}E_{i}$ are
already used, this means that $p$ separates $s_{k+1}$ from both $t_{k+1}$ and
the endpoints of the edges in $E_{k+1}$, a contradiction to $P_{k+1}$ being a
path in the solution. Hence $p$ uses an edge of $E_{k+1}$.

Since the sets $\bigcup_{1\leq i\leq k}P_{i}\smallsetminus\bigcup_{1\leq i\leq
k}E_{i}$ and $E_{k+1}$ have equal size, it follows that each curve of the
matching

$\bigcup_{1\leq i\leq k}P_{i}\smallsetminus\bigcup_{1\leq i\leq k}E_{i}$

uses precisely one edge of $E_{k+1}$. Since the endpoints of the matching are
fixed, they induce an order on the matching curves which determines precisely
which edge of $E_{k+1}$ is used by which curve.

Altogether, this shows that the solution to
$X_{k+1},(s_{1},t_{1}),\ldots,(s_{k+1},t_{k+1})$ is unique up to homeomorphism
and uses all edges $\bigcup_{1\leq i\leq k}E_{i}$.

∎

Figure 7: The number of edges around the terminals is crucial (cf. Remark 4).

######  Remark 4.

In a topological $\operatorname{DPP}$ instance, the number of edges around the
terminals is crucial. Even just relaxing the conditions on $X_{k}$ by having 2
edges instead of 1 edge around terminal $t_{2}$ allows a quite different
solution to the topological $\operatorname{DPP}$. This solution uses no edge
around $t_{k}$, one edge around each of $t_{3},t_{3},\ldots,t_{k-1}$, and the
two edges around $t_{2}$ (Figure 7 shows this for $k=4$).

######  Theorem 4.

Let $k\in\mathbb{N}\smallsetminus\\{0\\}$. The graph $G_{k}$ contains a vital
linkage.

###### Proof.

Let $P_{1},\ldots,P_{k}$ be the linkage from Remark 2. We argue that it is
vital. For $k=1$ and $k=2$, one can easily verify that $G_{k}$ has a unique
embedding. For $k\geq 2$, contracting an edge at $s_{1}$ suffices to make
$G_{k}$ $3$-connected. Since 3-connected planar graphs have unique embeddings
[14], the graph $G_{k}$ also has a unique embedding, and it suffices to
consider our previous embedding of $G_{k}$ (cf. Figure 2). Let $D$ be the
minimal disc containing the grid in $G_{k}$. The disc $D$ together with
$E(G_{k})$ is the space $X_{k}$. The paths $P_{1},\ldots,P_{k}$ thus give a
solution to the topological $\operatorname{DPP}$ instance
$X_{k},(s_{1},t_{1}),\ldots,(s_{k},t_{k})$, which by Lemma 1 is unique and
uses all edges in $E_{k}$. Thus any linkage
$P_{1}^{\prime},\ldots,P_{k}^{\prime}$ with the same pattern as
$P_{1},\ldots,P_{k}$ can differ from $P_{1},\ldots,P_{k}$ only inside the
grid. Thus for each $y\in[2^{k}-1]$ there is a subpath $Q_{y}^{\prime}$ of
some path of the solution $P_{1}^{\prime},\ldots,P_{k}^{\prime}$, such that
the endpoints of $Q_{y}^{\prime}$ are $(y^{\prime},1)$ and $(y,2^{k}-1)$ for
some $y^{\prime}\in[2^{k}-1]$. Hence the family
$(Q^{\prime}_{y})_{y\in[2^{k}-1]}$ is a linkage between the first column and
the last column of the grid.

Suppose that $P_{1}^{\prime},\ldots,P_{k}^{\prime}$ indeed differs from
$P_{1},\ldots,P_{k}$. Then at least one path $Q_{y}^{\prime}$ contains a
vertical edge $e$ in the grid. Hence the column of $e$ contains at most
$2^{k}-3$ vertices that are not used by $Q_{y}^{\prime}$ and, by Menger’s
Theorem [9], the remaining $2^{k}-2$ paths of the family cannot be routed, a
contradiction.

∎

_Proof of Theorem 3_ Theorem 3 immediately follows from Theorem 4 and Remark
1. ∎

## Acknowledgements

This research was partially supported by the Deutsche Forschungsgemeinschaft,
project Graphstrukturtheorie und algorithmische Anwendungen, AD 411/1-1 and
project Graphstrukturtheorie im Übersetzerbau, KR 4970/1-1. We thank Frédéric
Mazoit for valuable discussions, especially for inspiring Remark 4. We also
thank an anonymous reviewer for suggesting very elegant shortenings of our
construction and proof.

## References

  * [1] Isolde Adler, Stavros G. Kolliopoulos, Philipp K. Krause, Daniel Lokshtanov, Saket Saurabh, and Dimitrios M. Thilikos.  Tight Bounds for Linkages in Planar Graphs.  In Luca Aceto, Monika Henzinger, and Jiri Sgall, editors, ICALP (1), volume 6755 of Lecture Notes in Computer Science, pages 110–121. Springer, 2011.
  * [2] Isolde Adler, Stavros G. Kolliopoulos, Philipp Klaus Krause, Daniel Lokshtanov, Saket Saurabh, and Dimitrios M. Thilikos.  Irrelevant vertices for the planar Disjoint Paths Problem.  Journal of Combinatorial Theory, Series B, 122:815 – 843, 2017\.
  * [3] Jim Geelen, Tony Huynh, and Ronny B. Richter.  Explicit bounds for graph minors.  Journal of Combinatorial Theory, Series B, 132:80 – 106, 2018.
  * [4] Richard M. Karp.  On the Computational Complexity of Combinatorial Problems.  Networks, 5:45–68, 1975.
  * [5] Philipp K. Krause.  Graph Decomposition in Routing and Compilers.  PhD thesis, 2016.
  * [6] James F. Lynch.  The equivalence of theorem proving and the interconnection problem.  SIGDA Newsletter, 5(3):31–36, 1975.
  * [7] Jiří Matoušek, Eric Sedgwick, Martin Tancer, and Uli Wagner.  Untangling two systems of noncrossing curves.  Israel Journal of Mathematics, 212(1):37–79, 2016.
  * [8] Frédéric Mazoit.  A single exponential bound for the redundant vertex theorem on surfaces.  CoRR, abs/1309.7820, 2013.
  * [9] Karl Menger.  Zur allgemeinen Kurventheorie.  Fundamenta Mathematicae, 10(1):96–115, 1927.
  * [10] Bojan Mohar and Carsten Thomassen.  Graphs on Surfaces. Johns Hopkins series in the mathematical sciences. Johns Hopkins University Press, 2001.
  * [11] Neil Robertson and Paul D. Seymour.  Graph Minors. XIII. The Disjoint Paths Problem.  Journal of Combinatorial Theory, Series B, 63(1):65–110, 1995.
  * [12] Neil Robertson and Paul D. Seymour.  Graph Minors. XXI. Graphs with unique linkages.  Journal of Combinatorial Theory, Series B, 99(3):583–616, 2009\.
  * [13] Neil Robertson and Paul D. Seymour.  Graph Minors. XXII. Irrelevant vertices in linkage problems.  Journal of Combinatorial Theory, Series B, 102(2):530–563, 2012\.
  * [14] Hassler Whitney.  Congruent graphs and the connectivity of graphs.  American Journal of Mathematics, 54(1):150–168, 1932.


# A variational non-linear constrained model for the inversion of FDEM data

A Buccini1 and P Díaz de Alba2

###### Abstract

Reconstructing the structure of the soil using non-invasive techniques is a
very relevant problem in many scientific fields, like geophysics and
archaeology. This can be done, for instance, with the aid of Frequency Domain
Electromagnetic (FDEM) induction devices. Inverting FDEM data is a very
challenging inverse problem, as the problem is extremely ill-posed, i.e.,
sensible to the presence of noise in the measured data, and non-linear.
Regularization methods substitute the original ill-posed problem with a well-
posed one whose solution is an accurate approximation of the desired one. In
this paper we develop a regularization method to invert FDEM data. We propose
to determine the electrical conductivity of the ground by solving a
variational problem. The minimized functional is made up by the sum of two
term: the data fitting term ensures that the recovered solution fits the
measured data, while the regularization term enforces sparsity on the
Laplacian of the solution. The trade-off between the two terms is determined
by the regularization parameter. This is achieved by minimizing an
$\ell_{2}-\ell_{q}$ functional with $0<q\leq 2$. Since the functional we wish
to minimize is non-convex, we show that the variational problem admits a
solution. Moreover, we prove that, if the regularization parameter is tuned
accordingly to the amount of noise present in the data, this model induces a
regularization method. Some selected numerical examples on synthetic and real
data show the good performances of our proposal.

##  1 Introduction

In this paper we consider a severely ill-posed problem that arises in
geophysics, namely the reconstruction of the electrical conductivity of the
ground using a Frequency Domain Electromagnetic (FDEM) induction device. In
this work we consider a two-dimensional vertical section of the ground and
assume the magnetic permeability to be known and constant. The extension to
the three-dimensional case and unknown magnetic permeability is
straightforward and, for simplicity, we do not consider it here. The
considered problem is of the form

(1)
$\arg\min_{\Sigma\in\mathcal{X}}\left\|\mathcal{M}(\Sigma)-\mathcal{B}\right\|_{\mathcal{Y}}^{2},$

where $\mathcal{X}$ and $\mathcal{Y}$ are Hilbert spaces,
$\Sigma\in\mathcal{X}$ represents the electrical conductivity,
$\mathcal{B}\in\mathcal{Y}$ collects the (possibly noisy) data measured by the
FDEM device, $\left\|\cdot\right\|_{\mathcal{Y}}$ denotes the norm defined on
$\mathcal{Y}$, and $\mathcal{M}:\mathcal{X}\rightarrow\mathcal{Y}$ is a non-
linear function. We assume that (1) is ill-posed, i.e., the solution may not
be unique and its computation is very sensible to the presence of noise in the
measured data. Therefore, the naive solution of (1) is usually a poor
approximation of the exact electrical conductivity. To compute a meaningful
solution we need to resort to regularization methods; see, e.g., [18] and
references therein for a more detailed discussion.

We represent the two-dimensional ground as the semi-infinite rectangle
$[a,b]\times[0,\infty)$, where the first dimension is the horizontal space
dimension (intuitively the “sea-level”) and the second dimension represents
the depth. Our purpose is to reconstruct an “image” of the electrical
conductivity on this rectangle. This can be done using a FDEM induction device
called Ground Conductivity Meter (GCM). Its principle of operation is based on
an alternating electrical current which flows through a small electric wire
coil (the transmitter). A second coil (the receiver) is positioned at a fixed
distance from the first one, and the two coil axes may be aligned either
vertically or horizontally with respect to the subsurface. The transmitting
coil generates an electromagnetic (EM) field above the surface of the ground,
a portion of which propagates into it. This EM field, called the primary field
$H_{P}$, induces an alternating electrical current within the ground,
generating in turn a secondary EM field $H_{S}$, which propagates back to the
surface and the air above. The second wire coil acts as a receiver, measuring
the amplitude and phase components of the ratio between the primary and
secondary EM fields. The complex measurements obtained by a GCM depend on some
instrument settings, like the orientation of the dipoles, the frequency of the
alternating current, the inter-coil distance, and the height of the instrument
above the ground.

We provide below a brief discussion on the functional $\mathcal{M}$ presented
in (1). As we will see, the described model is one-dimensional, i.e., the
measured data depends exclusively on the electrical conductivity below the
measuring point.

We assume that $N$ equispaced measurement sets $\mathbf{b}_{j}$ are performed
in some interval $[a,b]$ (in our experiments we set $a=0$m and $b=10$m). Each
set is obtained by collecting several measurements using different
configurations of the FDEM device. Denoting a sampling of the electrical
conductivity of the ground by

$\mathbf{\Sigma}=[\boldsymbol{\sigma}_{1},\ldots,\boldsymbol{\sigma}_{N}],$

we can discretize the minimization problem (1) as follows

(2)
$\mathbf{\Sigma}=\arg\min_{\mathbf{\Sigma}}\sum_{j=1}^{N}\left\|M(\boldsymbol{\sigma}_{j})-\mathbf{b}_{j}\right\|_{2}^{2},$

where $\left\|\cdot\right\|_{2}$ is the Euclidean vector norm,
$\boldsymbol{\Sigma}=[\boldsymbol{\sigma}_{1},\ldots,\boldsymbol{\sigma}_{N}]\in{\mathbb{R}}^{n\times
N}$,
$\mathbf{B}=[\mathbf{b}_{1},\ldots,\mathbf{b}_{N}]\in{\mathbb{C}}^{m\times
N}$, and the vector function
$\mathbf{M}(\mathbf{\Sigma})=[M(\boldsymbol{\sigma}_{1}),\ldots,M(\boldsymbol{\sigma}_{N})]\in{\mathbb{C}}^{m\times
N}$ returns the readings predicted by the model in the same order they were
arranged in the vector $\mathbf{b}_{j}$. This can be done because the model is
one-dimensional, i.e., the values of the $j$th column of $\mathbf{B}$ depend
only on the entries of the $j$th column of $\mathbf{\Sigma}$.

Therefore, once discretized, the two-dimensional problem is reduced to $N$
independent one-dimensional problems.

Although this may seem beneficial, since this decoupling can be greatly useful
in lowering the computational complexity of the problem at hand, this is not
completely true. A natural approach for solving the two-dimensional problem
would be to solve each one-dimensional problem independently and then “stack”
the obtained solutions in order to get a two-dimensional image; see [14].
However, as we show in our numerical examples, in practice this approach may
lead to poor reconstructions. In fact, the problem is severely ill-posed and
in real-world applications the measured data is affected by noise, thus the
obtained solution can deviate (also substantially) from the exact
$\boldsymbol{\sigma}_{j}$ and this results in the stacked image to be
“spliced”. It is the purpose of this paper to propose a solution method that
couples the $\boldsymbol{\sigma}_{j}$ in order to avoid this “splicing” and to
counteract the ill-posedness of the problem. To this end, we consider the
following variational problem

(3) $\arg\min_{\mathbf{\Sigma}\geq
0}\frac{1}{2}\left\|\mathbf{M}(\mathbf{\Sigma})-\mathbf{B}\right\|_{F}^{2}+\frac{\gamma}{q}\left\|\mathbf{D}(\mathbf{\Sigma})\right\|_{q}^{q},$

where $\left\|\cdot\right\|_{F}$ denotes the Frobenius norm, $\gamma>0$, and
$0<q\leq 2$. With the notation $\mathbf{D}(\mathbf{\Sigma})$ we denote the
vector in ${\mathbb{R}}^{nN}$ defined by

$\mathbf{D}(\mathbf{\Sigma})=D(\mbox{{vec}}(\mathbf{\Sigma})),$

with ${D}={L}\otimes{I}+{I}\otimes{L}$, where $\otimes$ denotes the Kronecker
product, ${I}$ is the identity matrix, vec is the vectorization operator, and

${L}=\begin{bmatrix}2&-1\\\ -1&2&-1\\\ &\ddots&\ddots&\ddots\\\ &&-1&2&-1\\\
&&&-1&2\end{bmatrix}.$

The operator $\mathbf{D}$ (or equivalently the matrix $D$) is a discretization
of the Laplacian, thus we are considering a high order version of so-called
Total Variation regularization [29]. When $q=2$ the model (3) reduces to the
well-known Tikhonov minimization for non-linear problems, here with the
addition of the non-negativity constraint; see, e.g., [30] for a first
discussion on this kind of problems, [18, 22] for a recent discussion on
Tikhonov method and [1] for an algorithm for the solution of the non-
negatively constrained Tikhonov minimization for linear problems. When
$1<q<2$, the term $\left\|\mathbf{D}(\mathbf{\Sigma})\right\|_{q}^{q}$ is
convex and smooth, while for $q=1$ it is convex, but non-smooth. For $0<q<1$,
we define $\left\|\mathbf{x}\right\|_{q}^{q}=\sum_{j=1}^{n}|x_{j}|^{q}$ and we
still refer to this quantity as to $\ell_{q}$-norm, however, this is not a
norm since it does not satisfy the triangular inequality. Therefore, if
$0<q<1$, the second term in (3) is non-convex. Nevertheless, as shown in [25,
8, 5, 7, 9, 10], in imaging application it is beneficial to select $0<q<1$,
this is due to the fact that, in this case, the $\ell_{q}$-norm approximates
the $\ell_{0}$-norm, thus leading to reconstructions with sparse gradients.
This property is desirable in an approximate reconstruction.

Finally, observe that, since we have used zero boundary conditions for the
discretization of the Laplacian, it holds that $\mbox{{rank}}({D})=nN$.
Moreover, note that, since we know that the electrical conductivity is
positive, we imposed the non-negativity constraint.

The main contributions of this paper are the following. Firstly, we describe a
new variational model for computing an approximate solution of (1). The
novelty of the proposed model is the horizontal coupling of the columns of
$\mathbf{\Sigma}$ that allows us to obtain very accurate reconstructions of
the electrical conductivity of the ground that are not affected by the
“splicing” issue described above. To the best of our knowledge, this coupling
has never been considered before in the literature for this kind of problem.
Secondly, we prove the regularization properties of $\ell^{2}-\ell^{q}$
minimization for non-linear problems. This properties has been shown firstly
in [7] for linear problems and we extend this analysis to the non-linear case.
Finally, we propose an algorithm for the solution of the problem (3) that
extends to the non-linear case the one proposed in [25] that was constructed
for linear problems.

This paper is organized as follows: Section 2 describes the functional
$\mathcal{M}$ and its discretization, in Section 3 we show the theoretical
properties of (3). Section 4 provides an algorithm for the computation of an
approximate solution of (3), in Section 5 we show some numerical examples to
demonstrate the performances of the proposed approach, and we draw some
conclusions and outline future research in Section 6.

##  2 Modelization of the problem

We now briefly describe how the FDEM device can be modeled. The main non-
linear forward model which describes the interaction between the soil and the
FDEM induction device when the electrical conductivity and the magentic
permeability are known, has been described in [23].

In this model, the soil is assumed to have $n$ layers which are charaterized
by an electrical conductivity $\sigma_{i}$ (measured in S/m) and a magnetic
permeability $\mu_{i}$ (measured in H/m), for $i=1,\ldots,n$; see [15, 12] for
more details. The thickness of each layer, measured in meters, is denoted by
$d_{i}$, considering infinite the thikness of the deepest layer $d_{n}$.
Finally, the distance between the coils is represented by $\rho$ and the
height at which the measurements are taken by $h$.

Let us now consider the propagation constant
$u_{i}(\lambda)=\sqrt{\lambda^{2}+{\mathrm{i}}\sigma_{i}\mu_{i}\omega}$, where
${\mathrm{i}}$ denotes the imaginary unit, ${\rm Re}(u_{i}(\lambda))\geq 0$,
and $\omega$ is the angular frequency of the instrument, that is, $2\pi$ times
the frequency in Hertz. The variable $\lambda$ is the variable of integration
which ranges from zero to infinity measuring the ratio between the depth in
meters below the ground surface and the inter-coil distance $\rho$.

The surface admittance $Y_{i}(\lambda)$ at the top of each layer verifies the
recursion

(4)
$Y_{i}(\lambda)=N_{i}(\lambda)\frac{Y_{i+1}(\lambda)+N_{i}(\lambda)\tanh(d_{i}u_{i}(\lambda))}{N_{i}(\lambda)+Y_{i+1}(\lambda)\tanh(d_{i}u_{i}(\lambda))},$

for $i=n-1,\ldots,1$, where
$N_{i}(\lambda)=u_{i}(\lambda)/({\mathrm{i}}\mu_{i}\omega)$ represents the
characteristic admittance at the $i$-th layer; see [31]. Note that this
formulation of $Y_{i}$ is not computationally stable and to solve practical
problems different formulations have to be considered. However, our purpose
here is to describe the model and we do not dwell further on the computational
issues of the forward model. We refer the interested reader to [11] for a
discussion on the numerical implementation of the forward model. At the last
layer, i.e., when $i=n$, the characteristic admittance and the surface
admittance coincide being $Y_{n}(\lambda)=N_{n}(\lambda)$ used to initialize
the recursion (4). We remark here, that both the characteristic and the
surface admittances are functions of the frequency $\omega$ via the functions
$u_{i}(\lambda)$.

Now, the ratio of the secondary to the primary field for the vertical
($\nu=0$) and horizontal ($\nu=1$) orientation of the coils are given by

(5)
$M_{\nu}(\boldsymbol{\sigma},\boldsymbol{\mu};h,\omega,\rho)=-\rho^{3-\nu}\mathcal{H}_{\nu}\left[\lambda^{1-\nu}e^{-2h\lambda}R_{\omega,0}(\lambda)\right](\rho),\qquad\nu=0,1,$

where $\boldsymbol{\sigma}=(\sigma_{1},\ldots,\sigma_{n})^{T}$,
$\boldsymbol{\mu}=(\mu_{1},\ldots,\mu_{n})^{T}$,
$N_{0}(\lambda)=\lambda/({\mathrm{i}}\mu_{0}\omega)$, $\mu_{0}=4\pi\cdot
10^{-7}$H/m is the magnetic permeability of free space. The reflection factor
$R_{\omega,0}(\lambda)$ is defined by

$R_{\omega,0}(\lambda)=\frac{N_{0}(\lambda)-Y_{1}(\lambda)}{N_{0}(\lambda)+Y_{1}(\lambda)},$

with $Y_{1}(\lambda)$ computed by the recursion (4), and

$\mathcal{H}_{\nu}[f](\rho)=\int_{0}^{\infty}f(\lambda)J_{\nu}(\rho\lambda)\lambda\,d\lambda,\qquad\nu=0,1,$

is the Hankel transform, where $J_{0},J_{1}$ are first kind Bessel functions
of order 0 and 1, respectively.

Note that functions in (5) are complex valued functions. The imaginary part or
quadrature component of the field ratio is usually interpreted as the apparent
conductivity of the soil, while the magnetic permeability is related to the
real or in-phase component.

Simultaneous measurements with different inter-coil distances or different
operating frequencies can be recorded by recent FDEM devices at different
heights. We denote by
$\boldsymbol{\rho}=(\rho_{1},\ldots,\rho_{m_{\rho}})^{T}$,
$\mathbf{h}=(h_{1},\ldots,h_{m_{h}})^{T}$, and
$\boldsymbol{\omega}=(\omega_{1},\ldots,\omega_{m_{\omega}})^{T}$, the vectors
containing the loop-loop distances, the heights, and the angular frequencies
at which the readings were taken. We consider the corresponding
$m=2m_{\rho}m_{h}m_{\omega}$ data points $b^{\nu}_{tls}$, where
$t=1,\ldots,m_{\rho}$, $l=1,\ldots,m_{h}$, $s=1,\ldots,m_{\omega}$, while
$\nu\in\\{0,1\\}$ represents the vertical and horizontal orientations of the
coils, respectively. The observations $b^{\nu}_{tls}$ are rearranged in a
vector $\mathbf{b}\in{\mathbb{C}}^{m}$.

In various papers this non-linear model has been studied for different device
configurations and different techniques were applied; see, e.g., [23, 15, 17,
12, 13, 14]. An algorithm for the regularized inversion of this model has been
implemented in a Matlab package which includes a graphical user interface in
[11].

We remark that, for small values of the conductivity of the soil, a linear
model has been introduced in [26] and has been solved first in [3] and later
in [16] from the theoretical point of view, where an optimized solution method
has been proposed.

In the following, it is assumed that the contribution of the permeability
distribution to the overall response is negligible, i.e. $\mu=\mu_{0}$, so
that the measurements are considered to be sensitive merely to conductivity
values. However, in principle, the regularization approach discussed here can
be easily extended to include also the inversion for the $\boldsymbol{\mu}$
components; see [12].

##  3 Regularization property

In this section we consider the model (3) in the presence of noise, that is

(6) $\arg\min_{\mathbf{\Sigma}\geq
0}\frac{1}{2}\left\|\mathbf{M}(\mathbf{\Sigma})-\mathbf{B}^{\delta}\right\|_{F}^{2}+\frac{\gamma}{q}\left\|\mathbf{D}(\mathbf{\Sigma})\right\|_{q}^{q},$

where $\delta$ denotes the noise present in the data, i.e.,

$\left\|\mathbf{B}-\mathbf{B}^{\delta}\right\|\leq\delta.$

We would like to show that the model (6) induces a regularization method.
Namely, if $\gamma$ is chosen depending on the amount of noise that corrupts
the data, then the solutions of (6) converge to a minimum norm solution of the
noise free problem (3) as $\mathbf{B}^{\delta}\rightarrow\mathbf{B}$ or,
equivalently, $\delta\rightarrow 0$. Denote by $\mathcal{J}$ the functional
minimized in (6), i.e.,

(7)
$\mathcal{J}(\mathbf{\Sigma})=\frac{1}{2}\left\|\mathbf{M}(\mathbf{\Sigma})-\mathbf{B}^{\delta}\right\|_{F}^{2}+\frac{\gamma}{q}\left\|\mathbf{D}(\mathbf{\Sigma})\right\|_{q}^{q}+\iota_{0}(\mathbf{\Sigma}),$

where $\iota_{0}(\mathbf{\Sigma})$ is the indicator function of the non-
negative cone, namely

$\iota_{0}(\mathbf{\Sigma})=\left\\{\begin{array}[]{ll}0&\mbox{if
}\mathbf{\Sigma}\geq 0,\\\ \infty&\mbox{else.}\end{array}\right.$

Before showing the regularization property we need to first prove some
auxiliary results. The analysis of the model can be derived from the ones in
[7, 24], however, since some proofs are different, we report it here for the
convenience of the reader.

We first recall the following result.

######  Lemma 1 ([7]).

Let $\\{x_{j}\\}_{j\in{\mathbb{N}}}$ be a sequence of elements of
${\mathbb{R}}^{n}$ and let $q>0$. If the
$\left\|\mathbf{x}_{j}\right\|_{q}^{q}$ are uniformly bounded, i.e., if there
exists a constant $c>0$ independent of $j$ such that

$\left\|\mathbf{x}_{j}\right\|_{q}^{q}\leq c\quad\forall j\in{\mathbb{N}},$

then $\left\|\mathbf{x}_{j}\right\|_{2}^{2}$ is uniformly bounded.

We can now show our first preliminary result.

######  Proposition 2.

Let $\mathcal{J}$ be defined in (7), then $\mathcal{J}$ admits a global
minimizer.

######  Proof 3.

By definition $\mathcal{J}$ is lower semi-continuous and proper. Moreover,
since $\mbox{{rank}}(D)=nN$, it is easy to see that $\mathcal{J}$ is coercive.
In fact, assume that the sequence
$\\{\mathbf{\Sigma}_{j}\\}_{j\in{\mathbb{N}}}$ is such that
$\left\|\mathbf{\Sigma}_{j}\right\|_{F}^{2}\rightarrow\infty$ as
$j\rightarrow\infty$, then, since $\mbox{{rank}}(D)=nN$,
$\left\|\mathbf{D}(\mathbf{\Sigma}_{j})\right\|_{2}^{2}\rightarrow\infty$ as
$j\rightarrow\infty$. Consequently,
$\left\|\mathbf{D}(\mathbf{\Sigma}_{j})\right\|_{q}^{q}\rightarrow\infty$ as
$j\rightarrow\infty$ and thus
$\mathcal{J}(\mathbf{\Sigma}_{j})\rightarrow\infty$ as $j\rightarrow\infty$.

Since $\mathcal{J}$ is proper, there exists $\mathbf{\Sigma}$ such that
$\mathcal{J}(\mathbf{\Sigma})<\infty$ and we can define

$\varphi=\inf_{\mathbf{\Sigma}\in{\mathbb{R}}^{n\times
N}}\mathcal{J}(\mathbf{\Sigma}).$

By lower semi-continuity of $\mathcal{J}$, there exists a sequence
$\\{\mathbf{\Sigma}_{j}\\}_{j\in{\mathbb{N}}}$ and $M\geq 0$ such that

(8)
$\mathcal{J}(\mathbf{\Sigma}_{j})\rightarrow\varphi\;\mbox{as}\;j\rightarrow\infty\quad\mbox{and}\quad\mathcal{J}(\mathbf{\Sigma}_{j})\leq
M\;\forall j\in{\mathbb{N}}.$

In particular, $\left\|\mathbf{D}(\mathbf{\Sigma}_{j})\right\|_{q}^{q}\leq M$
for all $j\in{\mathbb{N}}$ and, thanks to Lemma 1, there exists
$\widetilde{M}$ such that
$\left\|\mathbf{D}(\mathbf{\Sigma}_{j})\right\|_{2}^{2}\leq\widetilde{M}$ for
all $j\in{\mathbb{N}}$. Since $\mbox{{rank}}({D})=nN$, the sequence
$\\{\mathbf{\Sigma}_{j}\\}_{j\in{\mathbb{N}}}$ is uniformly bounded and, thus,
it admits a convergent subsequence denoted by
$\\{\mathbf{\Sigma}_{j_{k}}\\}_{j_{k}\in{\mathbb{N}}}$, with
$\mathbf{\Sigma}_{j_{k}}\rightarrow\mathbf{\Sigma}^{*}$ as
$j_{k}\rightarrow\infty$. We would like to show that $\mathbf{\Sigma}^{*}$ is
a minimizer of $\mathcal{J}$. The definition of $\varphi$ yields

$\varphi\leq\mathcal{J}(\mathbf{\Sigma}^{*})\leq\liminf_{j_{k}\rightarrow\infty}\mathcal{J}(\mathbf{\Sigma}_{j_{k}})=\lim_{j_{k}\rightarrow\infty}\mathcal{J}(\mathbf{\Sigma}_{j_{k}})=\varphi,$

where the second to last equality follows from the lower semi-continuity of
$\mathcal{J}$ and the last one follows from (8). This shows that
$\mathcal{J}(\mathbf{\Sigma}^{*})=\varphi$, i.e., that $\mathbf{\Sigma}^{*}$
is a global minimizer of $\mathcal{J}$ which concludes the proof.

We are now in position to show our main result.

######  Theorem 1.

Let $0<q\leq 2$ be fixed and
$\\{\mathbf{B}^{\delta_{j}}\\}_{j\in{\mathbb{N}}}$ be a sequence such that,
for all $j\in{\mathbb{N}}$,
$\left\|\mathbf{B}-\mathbf{B}^{\delta_{j}}\right\|_{F}\leq\delta_{j}$. Assume
that $\delta_{j}\rightarrow 0$ as $j\rightarrow\infty$. Let
$\\{\gamma_{j}\\}_{j\in{\mathbb{N}}}$ be a sequence of positive real numbers
such that

$\gamma_{j}\rightarrow
0\quad\mbox{and}\quad\frac{\delta^{2}_{j}}{\gamma_{j}}\rightarrow
0\quad\mbox{as}\quad j\rightarrow\infty.$

Denote by $\mathcal{J}_{j}$ the functional

$\mathcal{J}_{j}(\mathbf{\Sigma})=\frac{1}{2}\left\|\mathbf{M}(\mathbf{\Sigma})-\mathbf{B}^{\delta_{j}}\right\|_{F}^{2}+\frac{\gamma_{j}}{q}\left\|\mathbf{D}(\mathbf{\Sigma})\right\|_{q}^{q}+\iota_{0}(\mathbf{\Sigma})$

and, for all $j\in{\mathbb{N}}$, let

$\mathbf{\Sigma}_{j}\in\arg\min_{\mathbf{\Sigma}}\mathcal{J}_{j}(\mathbf{\Sigma}).$

Then there exists a convergent subsequence of
$\\{\mathbf{\Sigma}_{j}\\}_{j\in{\mathbb{N}}}$, denoted by
$\\{\mathbf{\Sigma}_{j_{k}}\\}_{j_{k}\in{\mathbb{N}}}$, such that

$\mathbf{\Sigma}_{j_{k}}\rightarrow\mathbf{\Sigma}^{*}\quad\mbox{as}\quad
j_{k}\rightarrow\infty$

and

$\mathbf{\Sigma}^{*}\in\arg\min\\{\left\|\mathbf{D}(\mathbf{\Sigma})\right\|_{q}^{q}:\;\mathbf{M}(\mathbf{\Sigma})=\mathbf{B},\;\mathbf{\Sigma}\geq
0\\},$

assuming that this set is not empty.

######  Proof 4.

First let us observe that, thanks to Proposition 2, the sequence
$\\{\mathbf{\Sigma}_{j}\\}_{j\in{\mathbb{N}}}$ is well defined.

Since $\mathbf{\Sigma}_{j}$ is a global minimizer of $\mathcal{J}_{j}$, we
have that for all $\mathbf{\Sigma}$ it holds

$\mathcal{J}_{j}(\mathbf{\Sigma}_{j})\leq\mathcal{J}_{j}(\mathbf{\Sigma}).$

In particular, let
$\mathbf{\Sigma}^{\dagger}\in\arg\min\\{\left\|\mathbf{D}(\mathbf{\Sigma})\right\|_{q}^{q}:\;\mathbf{M}(\mathbf{\Sigma})=\mathbf{B}\\}$,
then

$\mathcal{J}_{j}(\mathbf{\Sigma}_{j})\leq\mathcal{J}_{j}(\mathbf{\Sigma}^{\dagger}).$

Recalling that $\mathbf{M}(\mathbf{\Sigma}^{\dagger})=\mathbf{B}$ and that
$\left\|\mathbf{B}-\mathbf{B}^{\delta_{j}}\right\|_{F}\leq\delta_{j}$, for all
$j\in{\mathbb{N}}$, we have that there exists $j_{0}$ such that, for all
$j>j_{0}$, it holds

$\displaystyle\frac{1}{2}\left\|\mathbf{M}(\mathbf{\Sigma}_{j})-\mathbf{B}^{\delta_{j}}\right\|_{F}^{2}+\frac{\gamma_{j}}{q}\left\|\mathbf{D}(\mathbf{\Sigma}_{j})\right\|_{q}^{q}$
$\displaystyle\leq\frac{1}{2}\left\|\mathbf{M}(\mathbf{\Sigma}^{\dagger})-\mathbf{B}^{\delta_{j}}\right\|_{F}^{2}+\frac{\gamma_{j}}{q}\left\|\mathbf{D}(\mathbf{\Sigma}^{\dagger})\right\|_{q}^{q}$
$\displaystyle\leq\frac{1}{2}\delta_{j}^{2}+\frac{\gamma_{j}}{q}\left\|\mathbf{D}(\mathbf{\Sigma}^{\dagger})\right\|_{q}^{q}$
$\displaystyle\leq C,$

where $C$ is a constant independent of $j$, where the last inequality follows
from the fact that $\delta_{j},\gamma_{j}\rightarrow 0$ as
$j\rightarrow\infty$, and we observed that
$\iota_{0}(\mathbf{\Sigma}_{j})=\iota_{0}(\mathbf{\Sigma}^{\dagger})=0$. In
particular, we have that, for all $j>j_{0}$,
$\left\|\mathbf{D}(\mathbf{\Sigma}_{j})\right\|_{q}^{q}\leq C$. Thanks to
$\mbox{{rank}}({D})=nN$ and Lemma 1, we have that
$\left\|\mathbf{\Sigma}_{j}\right\|_{F}\leq\widetilde{C}$ for a certain
constant $\widetilde{C}$. Since the sequence
$\\{\mathbf{\Sigma}_{j}\\}_{j\in{\mathbb{N}}}$ is uniformly bounded (for
$j>j_{0}$), it admits a converging subsequence
$\\{\mathbf{\Sigma}_{j_{k}}\\}_{j_{k}\in{\mathbb{N}}}$ and let
$\mathbf{\Sigma}^{*}$ denote its limit. We first show that
$\mathbf{M}(\mathbf{\Sigma}^{*})=\mathbf{B}$.

$\displaystyle 0$
$\displaystyle\leq\frac{1}{2}\left\|\mathbf{M}(\mathbf{\Sigma}^{*})-\mathbf{B}\right\|^{2}_{F}\leq\liminf_{j_{k}\rightarrow\infty}\frac{1}{2}\left\|\mathbf{M}(\mathbf{\Sigma}_{j_{k}})-\mathbf{B}^{\delta_{j_{k}}}\right\|^{2}_{F}$
$\displaystyle\leq\liminf_{j_{k}\rightarrow\infty}\frac{1}{2}\left\|\mathbf{M}(\mathbf{\Sigma}_{j_{k}})-\mathbf{B}^{\delta_{j_{k}}}\right\|^{2}_{F}+\frac{\gamma_{j_{k}}}{q}\left\|\mathbf{D}(\mathbf{\Sigma}_{j_{k}})\right\|_{q}^{q}+\iota_{0}(\mathbf{\Sigma}_{j_{k}})$
$\displaystyle\leq\liminf_{j_{k}\rightarrow\infty}\frac{\delta_{j_{k}}}{2}+\frac{\gamma_{j_{k}}}{q}\left\|\mathbf{D}(\mathbf{\Sigma}^{\dagger})\right\|_{q}^{q}+\iota_{0}(\mathbf{\Sigma}^{\dagger})=0,$

i.e., $\mathbf{M}(\mathbf{\Sigma}^{*})=\mathbf{B}$. We now show that
$\mathbf{\Sigma}^{*}$ minimizes
$\left\|\mathbf{D}(\mathbf{\Sigma}^{*})\right\|_{q}^{q}$. Recall that
$\left\|\mathbf{D}(\mathbf{\Sigma}^{\dagger})\right\|_{q}^{q}$ is minimum,
thus

$\displaystyle\frac{1}{q}\left\|\mathbf{D}(\mathbf{\Sigma}^{\dagger})\right\|_{q}^{q}$
$\displaystyle\leq\frac{1}{q}\left\|\mathbf{D}(\mathbf{\Sigma}^{*})\right\|_{q}^{q}\leq\liminf_{j_{k}\rightarrow\infty}\frac{1}{q}\left\|\mathbf{D}(\mathbf{\Sigma}_{j_{k}})\right\|_{q}^{q}+\iota_{0}(\mathbf{\Sigma}_{j_{k}})$
$\displaystyle\leq\liminf_{j_{k}\rightarrow\infty}\frac{1}{2\gamma_{j_{k}}}\left\|\mathbf{M}(\mathbf{\Sigma}_{j_{k}})-\mathbf{B}^{\delta_{j_{k}}}\right\|_{F}^{2}+\frac{1}{q}\left\|\mathbf{D}(\mathbf{\Sigma}_{j_{k}})\right\|_{q}^{q}+\iota_{0}(\mathbf{\Sigma}_{j_{k}})$
$\displaystyle=\liminf_{j_{k}\rightarrow\infty}\frac{1}{\gamma_{j_{k}}}\mathcal{J}_{j_{k}}(\mathbf{\Sigma}_{j_{k}})$
$\displaystyle\leq\liminf_{j_{k}\rightarrow\infty}\frac{1}{\gamma_{j_{k}}}\mathcal{J}_{j_{k}}(\mathbf{\Sigma}^{\dagger})$
$\displaystyle=\liminf_{j_{k}\rightarrow\infty}\frac{1}{2\gamma_{j_{k}}}\left\|\mathbf{M}(\mathbf{\Sigma}^{\dagger})-\mathbf{B}^{\delta_{j_{k}}}\right\|_{F}^{2}+\frac{1}{q}\left\|\mathbf{D}(\mathbf{\Sigma}^{\dagger})\right\|_{q}^{q}+\iota_{0}(\mathbf{\Sigma}^{\dagger})$
$\displaystyle\leq\liminf_{j_{k}\rightarrow\infty}\frac{\delta_{j_{k}}^{2}}{2\gamma_{j_{k}}}+\frac{1}{q}\left\|\mathbf{D}(\mathbf{\Sigma}^{\dagger})\right\|_{q}^{q}+\iota_{0}(\mathbf{\Sigma}^{\dagger})=\frac{1}{q}\left\|\mathbf{D}(\mathbf{\Sigma}^{\dagger})\right\|_{q}^{q}.$

We have shown that
$\left\|\mathbf{D}(\mathbf{\Sigma}^{*})\right\|_{q}^{q}=\left\|\mathbf{D}(\mathbf{\Sigma}^{\dagger})\right\|_{q}^{q}$.
Moreover, since the non-negative cone is a closed set it is obvious that
$\mathbf{\Sigma}^{*}\geq 0$, this coupled with the fact that
$\mathbf{M}(\mathbf{\Sigma}^{*})=\mathbf{B}$, concludes the proof.

##  4 Minimizing algorithm

We wish now to discuss how to compute a solution of (6). Consider the
minimization problem

(9)
$\arg\min_{\mathbf{\Sigma},\mathbf{\Xi}}\frac{1}{2}\left\|\mathbf{M}(\mathbf{\Sigma})-\mathbf{B}^{\delta}\right\|_{F}^{2}+\frac{\gamma}{q}\left\|\mathbf{D}(\mathbf{\Xi})\right\|_{q}^{q}+\frac{\beta}{2}\left\|\mathbf{\Sigma}-\mathbf{\Xi}\right\|_{F}^{2}+\iota_{0}(\mathbf{\Sigma}),$

where $\beta>0$ is a fixed parameter. Obviously, if $\beta$ is large enough,
the solutions of (6) and (9) are the same. This reformulation is commonly
performed in optimization so that the alternating minimization algorithm can
be used. By applying alternating minimization to solve (9), we obtain the
following iterations

(10)
$\left\\{\begin{array}[]{l}\displaystyle\mathbf{\Sigma}^{(k+1)}=\arg\min_{\mathbf{\Sigma}}\frac{1}{2}\left\|\mathbf{M}(\mathbf{\Sigma})-\mathbf{B}^{\delta}\right\|_{F}^{2}+\frac{\beta}{2}\left\|\mathbf{\Sigma}-\mathbf{\Xi}^{(k)}\right\|_{F}^{2}+\iota_{0}(\mathbf{\Sigma})\\\
\displaystyle\mathbf{\Xi}^{(k+1)}=\arg\min_{\mathbf{\Xi}}\frac{\gamma}{q}\left\|\mathbf{D}(\mathbf{\Xi})\right\|_{q}^{q}+\frac{\beta}{2}\left\|\mathbf{\Sigma}^{(k+1)}-\mathbf{\Xi}\right\|_{F}^{2}.\end{array}\right.$

The convergence of the iterations in (10) is guaranteed by the results in
[21]. In particular, it holds the following

######  Theorem 2.

Let $\mathbf{\Sigma}^{(k)}$ and $\mathbf{\Xi}^{(k)}$ denote the iterates
defined in (10), then

$(\mathbf{\Sigma}^{(k)},\mathbf{\Xi}^{(k)})\rightarrow(\mathbf{\Sigma}^{*},\mathbf{\Xi}^{*})\quad\mbox{as}\quad
k\rightarrow\infty,$

where $(\mathbf{\Sigma}^{*},\mathbf{\Xi}^{*})$ is a stationary point of the
functional in (9) that depends on the initial guess $\mathbf{\Xi}^{(0)}$.
Moreover, if $\beta$ is large enough (possibly $\beta=\infty$),
$\mathbf{\Sigma}^{*}=\mathbf{\Xi}^{*}=\widehat{\mathbf{\Sigma}}$, where
$\widehat{\mathbf{\Sigma}}$ is a stationary point of the functional in (6)
that depends on $\mathbf{\Xi}^{(0)}$.

Remark.Note that, since the minimized functional in (9) is non-convex it might
have multiple global and local minima. Moreover, it may have several saddle
points. However, in general, determining the global minimum of a non-convex
function is a NP-hard problem. Therefore, in optimization theory, is usually
considered a good enough result determining a stationary point of the
minimized functional. Indeed, our numerical results show that the proposed
approach is able to provide accurate reconstructions.

###  4.1 Implementation details

In the following, we detail how to numerically solve the two minimization
subproblems in (10) at each iteration.

For the solution of the $\mathbf{\Sigma}$ subproblem in (10) let us first
observe that we can write

$\mathbf{M}(\mathbf{\Sigma})=[M(\boldsymbol{\sigma}_{1}),\ldots,M(\boldsymbol{\sigma}_{N})],$

where $M(\boldsymbol{\sigma}_{j})$ is a column vector. We can then rewrite the
subproblem as

$\mathbf{\Sigma}^{(k+1)}=\arg\min_{\mathbf{\Sigma}}\sum_{j=1}^{N}\left[\frac{1}{2}\left\|M(\boldsymbol{\sigma}_{j})-\mathbf{b}_{j}^{\delta}\right\|_{2}^{2}+\frac{\beta}{2}\left\|\boldsymbol{\sigma}_{j}-\boldsymbol{\xi}^{(k)}_{j}\right\|_{2}^{2}+\iota_{0}(\boldsymbol{\sigma}_{j})\right],$

where $\mathbf{b}_{j}^{\delta}$ and $\boldsymbol{\xi}^{(k)}_{j}$ denote the
$j$th column of $\mathbf{B}^{\delta}$ and $\mathbf{\Xi}^{(k)}$, respectively,
and $\boldsymbol{\sigma}_{j}$ is the $j$th column of $\mathbf{\Sigma}$. Thus,
we can write

$\boldsymbol{\sigma}_{j}^{(k+1)}=\arg\min_{\boldsymbol{\sigma}}\frac{1}{2}\left\|M(\boldsymbol{\sigma})-\mathbf{b}_{j}^{\delta}\right\|_{2}^{2}+\frac{\beta}{2}\left\|\boldsymbol{\sigma}-\boldsymbol{\xi}^{(k)}_{j}\right\|_{2}^{2}+\iota_{0}(\boldsymbol{\sigma}),\quad
j=1,\ldots,N,$

i.e., the $\mathbf{\Sigma}$ subproblem decouples in $N$ independent one-
dimensional subproblems. This problem can be rewritten as

(11)
$\boldsymbol{\sigma}_{j}^{(k+1)}=\arg\min_{\boldsymbol{\sigma}}\frac{1}{2}\left\|\begin{bmatrix}M(\boldsymbol{\sigma})\\\
\sqrt{\beta}\boldsymbol{\sigma}\end{bmatrix}-\begin{bmatrix}\mathbf{b}_{j}^{\delta}\\\
\sqrt{\beta}\boldsymbol{\xi}^{(k)}_{j}\end{bmatrix}\right\|_{2}^{2}+\iota_{0}(\boldsymbol{\sigma}),\quad
j=1,\ldots,N.$

Observe that every $\boldsymbol{\sigma}_{j}^{(k+1)}$ subproblems are
independent and thus can be solved in parallel. For the solution of the
$\boldsymbol{\sigma}_{j}^{(k+1)}$, we consider a slight modification of the
algorithm proposed in [15, 17, 12, 13, 14, 11]. Here we outline the algorithm
and describe the modifications we made.

Let us first consider the following non-linear least-squares problem for a
single column $M(\boldsymbol{\sigma}_{j})$ and $\mathbf{b}_{j}^{\delta}$ of
$\mathbf{M}(\mathbf{\Sigma})$ and $\mathbf{B}^{\delta}$, respectively, i.e.

(12)
$\min_{\boldsymbol{\sigma}}\frac{1}{2}\left\|{M}(\boldsymbol{\sigma})-\mathbf{b}_{j}^{\delta}\right\|_{2}^{2}+\iota_{0}(\boldsymbol{\sigma}).$

We solve problem (12) by the Guass–Newton method.

We denote by
$\mathbf{r}(\boldsymbol{\sigma})={M}(\boldsymbol{\sigma})-\mathbf{b}^{\delta}_{j}$
the complex residual vector as a function of the conductivity
$\boldsymbol{\sigma}$. At each step of the iterative algorithm we minimize the
2-norm of a linear approximation of the residual and compute
$\boldsymbol{\sigma}^{(l+1)}=\boldsymbol{\sigma}^{(l)}+\mathbf{q}^{(l)}$,
where

(13)
$\displaystyle\mathbf{q}^{(l)}=\arg\min_{\mathbf{q}}\left\|\mathbf{r}(\boldsymbol{\sigma}^{(l)})+J(\boldsymbol{\sigma}^{(l)})\mathbf{q}\right\|_{2},$

and $J(\boldsymbol{\sigma}^{(l)})$ is the Jacobian matrix of the function
$\mathbf{r}$ computed in $\boldsymbol{\sigma}^{(l)}$.

Being the residual function $\mathbf{r}$ complex-valued, we solve problem (13)
by stacking the real and imaginary part of the residual as follows (see [14])

$\widetilde{\mathbf{r}}(\boldsymbol{\sigma})=\begin{bmatrix}\mbox{Re}(\mathbf{r}(\boldsymbol{\sigma}))\\\
\mbox{Im}(\mathbf{r}(\boldsymbol{\sigma}))\end{bmatrix}\in{\mathbb{R}}^{2m},\quad\widetilde{J}(\boldsymbol{\sigma})=\begin{bmatrix}\mbox{Re}(J(\boldsymbol{\sigma}))\\\
\mbox{Im}(J(\boldsymbol{\sigma}))\end{bmatrix}\in{\mathbb{R}}^{2m\times n}.$

In the same way we set $\widetilde{\mathbf{r}}$ and $\widetilde{J}$, we
rearrange the vectors $\mathbf{M}$ and $\mathbf{b}^{\delta}$, and we denote
them by $\mathbf{\widetilde{M}}$ and $\widetilde{\mathbf{b}^{\delta}}$,
respectively. So, we replace (13) by

(14)
$\min_{\mathbf{q}}\left\|\widetilde{\mathbf{r}}(\boldsymbol{\sigma}^{(l)})+\widetilde{J}_{l}\mathbf{q}\right\|_{2},$

with $\widetilde{J}_{l}=\widetilde{J}(\boldsymbol{\sigma}^{(l)})$, and the
iterative method becomes

$\boldsymbol{\sigma}^{(l+1)}=\boldsymbol{\sigma}^{(l)}+\alpha_{l}\mathbf{q}^{(l)}=\boldsymbol{\sigma}^{(l)}-\alpha_{l}\widetilde{J}_{l}^{\dagger}\,\widetilde{\mathbf{r}}(\boldsymbol{\sigma}^{(l)}),$

where $\widetilde{J}_{l}^{\dagger}$ is the Moore–Penrose pseudoinverse of
$\widetilde{J}_{l}$ and $\alpha_{l}$ is a damping parameter which ensures the
convergence. This parameter is determined by coupling the Armijo–Goldstein
principle [2] to the positivity constraint $\boldsymbol{\sigma}^{(l+1)}\geq 0$
(see [12, 15]) so that

(15)
$\left\|\mathbf{\widetilde{r}}(\mathbf{\sigma}^{(l)})\right\|_{2}^{2}-\left\|\mathbf{\widetilde{r}}(\mathbf{\sigma}^{(l)}+\alpha_{l}\mathbf{q}^{(l)})\right\|_{2}^{2}\geq\frac{1}{2}\alpha_{l}\left\|\widetilde{J}_{l}\mathbf{q}^{(l)}\right\|_{2}^{2}\quad\mbox{with}\quad\boldsymbol{\sigma}^{(l+1)}\geq
0$

is verified. The analytical expression of the Jacobian matrices with respect
to the electrical conductivity and the magnetic permeability were computed in
[15] and [12], respectively, where it has also been proved that the
computation of the analitical expression of the Jacobian matrices is faster
than their finite difference approximations.

It is well known that the minimization problem (14) is extremely ill-posed,
meaning that the matrix $\widetilde{J}_{l}$ is severely ill-conditioned. In
order to overcome this difficulty, we apply the generalized truncated singular
value decomposition (GTSVD), as in [15, 17, 12, 14, 11], for stably computing
an approximation of $\widetilde{J}_{l}^{\dagger}$. We introduce a
regularization matrix $\widehat{L}\in{\mathbb{R}}^{p\times 2n}$ ($p\leq 2n$),
whose null space approximately contains the sought solution [28]. Under the
assumption
$\mathcal{N}(\widetilde{J}_{l})\cap\mathcal{N}(\widehat{L})=\\{\mathbf{0}\\}$,
problem (14) is replaced by

(16)
$\min_{\mathbf{q}\in\mathcal{S}}\left\|L\mathbf{q}\right\|_{2}^{2},\qquad\mathcal{S}=\\{\mathbf{q}\in{\mathbb{R}}^{2n}\,:\,\widetilde{J}_{l}^{T}\widetilde{J}_{l}\mathbf{q}=-\widetilde{J}_{l}^{T}\widetilde{\mathbf{r}}(\boldsymbol{\sigma}^{(l)})\\}.$

Note that, thanks to the assumption above on the null spaces of
$\widetilde{J}_{l}$ and $\widehat{L}$, the solution of (16) is unique. Very
common choices for $L$ are the discretization of the first or second
derivative operators.

Let the generalized singular value decomposition (GSVD) [20] of the matrix
pair $(\widetilde{J}_{l},\widehat{L})$ be

$\widetilde{J}_{l}=U\Sigma_{J}Z^{-1},\qquad\widehat{L}=V\Sigma_{L}Z^{-1},$

where $U$ and $V$ are matrices with orthonormal columns $\mathbf{u}_{i}$ and
$\mathbf{v}_{i}$, respectively, $Z$ is a non-singular matrix with columns
$\mathbf{z}_{i}$, and $\Sigma_{J}$, $\Sigma_{L}$ are diagonal matrices with
diagonal entries $c_{i}$ and $s_{i}$, which are the singular values of
$\widetilde{J}$ and $\widehat{L}$, respectively. Under the assumption that
$m=2m_{\rho}m_{h}m_{\omega}<2n$, common in the generality of cases, the
truncated GSVD (TGSVD) solution $\mathbf{q}^{(l)}$ (see [22] for details) can
be written as

(17) $\mathbf{q}^{(l)}_{\rm
GSVD}=-\sum_{i=p-\ell+1}^{p}\frac{\mathbf{u}_{i}^{T}\widetilde{\mathbf{r}}^{(l)}}{c_{i-2n+\kappa}}\,\mathbf{z}_{i}-\sum_{i=p+1}^{2n}(\mathbf{u}_{i}^{T}\widetilde{\mathbf{r}}^{(l)})\,\mathbf{z}_{i},$

where $\ell$ is fixed, $\kappa=\mathop{\operator@font
rank}\nolimits(\widetilde{J}_{l})$, $\ell=1,\ldots,\kappa+p-2n$ is the
regularization parameter, and
$\mathbf{r}^{(l)}=\mathbf{r}(\boldsymbol{\sigma}^{(l)})$.

The resulting regularized damped Gauss–Newton method reads

$\boldsymbol{\sigma}^{(l+1)}=\boldsymbol{\sigma}^{(l)}+\alpha_{l}\mathbf{q}^{(l)}_{\rm
GSVD},$

with $\ell$ fixed and $\alpha_{l}$ determined at each step as in (15).

We now apply this algorithm to our specific case. In particular, we have to
solve (11). We apply the previously described procedure to the following
function getting

(18)
$\widehat{M}(\boldsymbol{\sigma})=\begin{bmatrix}\widetilde{M}(\boldsymbol{\sigma})\\\
\sqrt{\beta}\boldsymbol{\sigma}\end{bmatrix},\quad\widehat{J}(\boldsymbol{\sigma})=\begin{bmatrix}\widetilde{J}(\boldsymbol{\sigma})\\\
\sqrt{\beta}I\end{bmatrix},$

where $I$ is the identity matrix. We remind here that $\widehat{M}$ is not a
complex valued function anymore since we already stacked the real and
imaginary parts of $M$ and $J$.

Let $\boldsymbol{\sigma}_{j}^{(k,l)}$ be an approximation of
$\boldsymbol{\sigma}_{j}^{(k+1)}$. Denote by $\widehat{J}_{l}$ the Jacobian of
$\widehat{M}$ in $\boldsymbol{\sigma}_{j}^{(k,l)}$ and consider the GSVD of
the pair $\\{\widehat{J}_{l},\widehat{L}\\}$

$\widehat{J}_{l}=\widehat{U}\widehat{\Sigma}_{J}\widehat{Z}^{-1},\qquad\widehat{L}=\widehat{V}\widehat{\Sigma}_{L}\widehat{Z}^{-1}.$

Let

$\widehat{\mathbf{r}}^{(l)}=\widehat{M}(\boldsymbol{\sigma}_{j}^{(k,l)})-\begin{bmatrix}\mathbf{b}_{j}^{\delta}\\\
\boldsymbol{\xi}^{(k)}_{j}\end{bmatrix}=\begin{bmatrix}M(\boldsymbol{\sigma}_{j}^{(k,l)})-\mathbf{b}_{j}^{\delta}\\\
\sqrt{\beta}(\boldsymbol{\sigma}_{j}^{(k,l)}-\boldsymbol{\xi}^{(k)}_{j})\end{bmatrix},$

then we compute, analogously to (17),

$\mathbf{q}^{(l)}_{\rm
GSVD}=-\sum_{i=p-\ell+1}^{p}\frac{\mathbf{u}_{i}^{T}\widehat{\mathbf{r}}^{(l)}}{c_{i-2n+\kappa}}\,\mathbf{z}_{i}-\sum_{i=p+1}^{2n}(\mathbf{u}_{i}^{T}\widehat{\mathbf{r}}^{(l)})\,\mathbf{z}_{i},$

leading to the iteration

$\boldsymbol{\sigma}_{j}^{(k+1,l)}=\boldsymbol{\sigma}_{j}^{(k,l)}+\alpha_{l}\mathbf{q}^{(l)}_{\rm
GSVD},$

where $\alpha_{l}$ is determined by the Armijo-Goldstein rule so that

$\left\|\widehat{\mathbf{r}}^{(l)}\right\|^{2}-\left\|\widehat{\mathbf{r}}^{(l+1)}\right\|_{2}^{2}\geq\frac{1}{2}\alpha_{l}\left\|\widehat{J}_{l}\mathbf{q}^{(l)}_{\rm
GSVD}\right\|_{2}^{2}\quad\mbox{and}\quad\boldsymbol{\sigma}^{(k,l+1)}\geq 0.$

We now move to the $\mathbf{\Xi}$ subproblem. To solve this problem we
consider the majorization-minimization algorithm proposed by Huang et al. in
[25] and furtherly developed in [5, 8, 9]. We briefly describe the algorithm
presented in [25]. To simplify the computations we consider a modified version
of the $\mathbf{D}$ operator, namely we impose reflexive boundary conditions
to the discretization of the Laplacian. In this way the obtained matrix has an
exploitable structure that helps in the computations. In detail, with abuse of
notation, we write

${L}=\begin{bmatrix}1&-1\\\ -1&2&-1\\\ &\ddots&\ddots&\ddots\\\ &&-1&2&-1\\\
&&&-1&1\end{bmatrix}.$

Thanks to the structure of $L$ the matrix $D$, which is defined by $D=L\otimes
I+I\otimes L$, is the sum of a block Hankel with Hankel blocks matrix, a block
Toeplitz with Hankel blocks matrix, a block Hankel with Toeplitz blocks
matrix, and a block Toeplitz with Toeplitz blocks matrix. We recall that
Toeplitz matrices are matrices that are constant on the diagonals and Hankel
matrices are matrices that are constant on the anti-diagonals. Since
$[-1\;2\;-1]$ is symmetric, if ${C}$ denotes the discrete cosine transform
matrix, we have that

(19) ${D}={C}^{T}{\Lambda}{C},$

where ${\Lambda}$ is a diagonal matrix. The diagonal elements of ${\Lambda}$
are computed as the cosine coefficients of the first column of ${D}$; see [27]
for more details.

We now describe the Majorization-Minimization (MM) procedure that we use for
solving the $\mathbf{\Xi}$ subproblem. Let us first rewrite the minimization
problem in (10) as follows

(20)
$\mbox{{vec}}(\mathbf{\Xi}^{(k+1)})=\arg\min_{\boldsymbol{\xi}}\frac{1}{2}\left\|\boldsymbol{\xi}-\boldsymbol{\sigma}^{(k+1)}\right\|_{2}^{2}+\frac{\gamma}{q\beta}\left\|{D}\boldsymbol{\xi}\right\|_{q}^{q},$

where $\boldsymbol{\sigma}^{(k+1)}=\mbox{{vec}}(\mathbf{\Sigma}^{(k+1)})$. The
MM algorithm generates a sequence of vectors $\boldsymbol{\xi}^{(k,l)}$ that
converges to an approximate solution of (20). Firstly, if $q\leq 1$, we need
to smooth the $q$-norm so that it is differentiable. Let $\varepsilon>0$ be a
small constant, then, for $\mathbf{x}\in{\mathbb{R}}^{n}$, it holds

$\left\|\mathbf{x}\right\|_{q}^{q}=\sum_{j=1}^{n}|x_{j}|^{q}\approx\sum_{j=1}^{n}\left(x_{j}^{2}+\varepsilon^{2}\right)^{q/2}=:\left\|\mathbf{x}\right\|_{q,\varepsilon}^{q}.$

Note that the function
$\mathbf{x}\mapsto\left\|\mathbf{x}\right\|_{q,\varepsilon}^{q}$ is
differentiable everywhere. We consider the smoothed problem

(21)
$\mbox{{vec}}(\mathbf{\Xi}^{(k+1)})=\arg\min_{\boldsymbol{\xi}}\frac{1}{2}\left\|\boldsymbol{\xi}-{\boldsymbol{\sigma}}^{(k+1)}\right\|_{2}^{2}+\frac{\gamma}{q\beta}\left\|{D}\boldsymbol{\xi}\right\|_{q,\varepsilon}^{q}=\arg\min_{\boldsymbol{\xi}}\widehat{\mathcal{J}}_{\varepsilon}(\mathbf{\boldsymbol{\xi}}).$

As pointed out in [8], the solutions of (20) and (21) are extremely similar
and this substitution does not have any negative effect.

Let $\boldsymbol{\xi}^{(k,l)}$ be an approximate solution of (21). We first
construct a quadratic tangent majorant of
$\widehat{\mathcal{J}}_{\varepsilon}$ in $\boldsymbol{\xi}^{(k,l)}$ that
majorizes it, i.e., a quadratic function
$\mathcal{Q}(\boldsymbol{\xi},\boldsymbol{\xi}^{(k,l)})$ such that

  * •

$\mathcal{Q}(\boldsymbol{\xi},\boldsymbol{\xi}^{(k,l)})\geq\widehat{\mathcal{J}}_{\varepsilon}(\boldsymbol{\xi})$
for all $\boldsymbol{\xi}$;

  * •

$\mathcal{Q}(\boldsymbol{\xi}^{(k,l)},\boldsymbol{\xi}^{(k,l)})=\widehat{\mathcal{J}}_{\varepsilon}(\boldsymbol{\xi}^{(k,l)})$;

  * •

$\nabla\mathcal{Q}(\boldsymbol{\xi}^{(k,l)},\boldsymbol{\xi}^{(k,l)})=\nabla\widehat{\mathcal{J}}_{\varepsilon}(\boldsymbol{\xi}^{(k,l)})$.

In [25] the authors provide two different choices for the construction of such
a functional, here we considered the so-called fixed approach. Denote by
$\widetilde{\mathbf{u}}^{(l)}=D\boldsymbol{\xi}^{(k,l)}$, then we compute the
vector $\mathbf{u}^{(l)}$ as

$\mathbf{u}^{(l)}=\widetilde{\mathbf{u}}^{(l)}\left(1-\left(\frac{(\widetilde{\mathbf{u}}^{(l)})^{2}+\varepsilon^{2}}{\varepsilon^{2}}\right)^{q/2-1}\right),$

where all the operations are meant element-wise. It is possible to see that
the function

$\mathcal{Q}(\boldsymbol{\xi},\boldsymbol{\xi}^{(k,l)})=\frac{1}{2}\left\|\boldsymbol{\xi}-\boldsymbol{\sigma}^{(k+1)}\right\|_{2}^{2}+\frac{\gamma\varepsilon^{q-2}}{2\beta}\left(\left\|{D}\boldsymbol{\xi}\right\|_{2}^{2}-2\langle\mathbf{u}^{(l)},{D}\boldsymbol{\xi}\rangle\right)+c,$

where $c$ is a constant independent of $\boldsymbol{\xi}$, is a quadratic
tangent majorant; see [25] for a derivation. An improved approximation of the
solution of (21) can be obtained as the unique minimizer of
$\mathcal{Q}(\boldsymbol{\xi},\boldsymbol{\xi}^{(k,l)})$, i.e.,

$\displaystyle\boldsymbol{\xi}^{(k,l+1)}$
$\displaystyle=\arg\min_{\boldsymbol{\xi}}\frac{1}{2}\left\|\boldsymbol{\xi}-{\boldsymbol{\sigma}}^{(k+1)}\right\|_{2}^{2}+\frac{\gamma\varepsilon^{q-2}}{2\beta}\left(\left\|D\boldsymbol{\xi}\right\|_{2}^{2}-2\langle\mathbf{u}^{(l)},{D}\boldsymbol{\xi}\rangle\right)$
$\displaystyle=\arg\min_{\boldsymbol{\xi}}\frac{1}{2}\left\|\begin{bmatrix}\mathbf{I}\\\
\sqrt{\eta}{D}\end{bmatrix}\boldsymbol{\xi}-\begin{bmatrix}{\boldsymbol{\sigma}^{(k+1)}}\\\
\sqrt{\eta}\mathbf{u}^{(l)}\end{bmatrix}\right\|_{2}^{2},$

where $\eta=\frac{\gamma\varepsilon^{q-2}}{\beta}$. Writing the normal
equation yields

$({I}+\eta{D}^{T}{D})\boldsymbol{\xi}^{(k,l+1)}=\boldsymbol{\sigma}^{(k+1)}+\eta
D^{T}\mathbf{u}^{(l)}.$

Using the factorization (19) we get

${C}^{T}\left({I}+\eta{\Lambda}^{T}{\Lambda}T\right){C}\boldsymbol{\xi}^{(k,l+1)}=\boldsymbol{\sigma}^{(k+1)}+\eta{D}^{T}\mathbf{u}^{(l)}.$

Thus, we can compute $\boldsymbol{\xi}^{(k,l+1)}$ as

$\boldsymbol{\xi}^{(k,l+1)}={C}^{T}\left({I}+\eta\Lambda^{T}\Lambda\right)^{-1}{C}(\boldsymbol{\sigma}^{(k+1)}+\eta{D}^{T}\mathbf{u}^{(l)}),$

where the inversion is well-defined since $\eta>0$. These computations can be
performed fairly inexpensively, since the system to solve is a diagonal one
and the application of the $n\times n$ cosine matrix can be performed in
$O(n\log n)$ operations by means of the dct algorithm.

We summarize all the computations in Algorithm 1.

input : $\mathbf{B}^{\delta}$, $\mathbf{\Xi}^{(0)}$, $\gamma>0$, $\beta>0$,
$0<q\leq 2$, $\ell\ll nN$, $\varepsilon>0$, $\sigma_{0}>0$.

1 $\eta=\frac{\gamma\varepsilon^{q-2}}{\beta}$;

2 for  _$k=0,1,2,\ldots$_ do

3 parfor  _$j=1,\ldots,N$_

4 $\boldsymbol{\sigma}^{(k,1)}=\sigma_{0}\mathbf{1}$;

5 for  _$l=1,2,\ldots$_ do

6
$\widehat{J}_{l}=\begin{bmatrix}{\widetilde{J}}(\boldsymbol{\sigma}^{(k,l)})\\\
\sqrt{\beta}I\end{bmatrix}$;

7 Compute the GSVD
${\widehat{J}}=\widehat{U}\widehat{\Sigma}_{J}\widehat{Z}^{-1},\;\widehat{L}=\widehat{V}\widehat{\Sigma}_{L}\widehat{Z}^{-1}$;

8
${\mathbf{r}}^{(l)}=\begin{bmatrix}M(\boldsymbol{\sigma}_{j}^{(k,l)})-\mathbf{b}_{j}^{\delta}\\\
\sqrt{\beta}(\boldsymbol{\sigma}_{j}^{(k,l)}-\boldsymbol{\xi}^{(k)}_{j})\end{bmatrix}$;

9 $\displaystyle\mathbf{q}^{(l)}_{\rm
GSVD}=-\sum_{i=p-\ell+1}^{p}\frac{\mathbf{u}_{i}^{T}\mathbf{r}^{(l)}}{c_{i-2n+\kappa}}\,\mathbf{z}_{i}-\sum_{i=p+1}^{2n}(\mathbf{u}_{i}^{T}\mathbf{r}^{(l)})\,\mathbf{z}_{i}$;

10 Find $\alpha_{l}$ s.t.
$\left\\{\begin{array}[]{l}\left\|\mathbf{r}^{(l)}\right\|^{2}-\left\|\begin{bmatrix}M\left(\mathbf{\sigma}^{(k,l)}+\alpha_{l}\mathbf{q}_{\rm
GSVD}^{(l)}\right)-\mathbf{b}_{j}^{\delta}\\\
\sqrt{\beta}(\boldsymbol{\sigma}^{(k,l)}+\alpha_{l}\mathbf{q}_{\rm
GSVD}^{(l)}-\boldsymbol{\xi}^{(k)}_{j})\end{bmatrix}\right\|^{2}\geq\frac{1}{2}\alpha_{l}\left\|\widehat{J}_{l}\mathbf{q}_{\rm
GSVD}^{(l)}\right\|^{2}\\\
\boldsymbol{\sigma}^{(k,l)}+\alpha_{l}\mathbf{q}_{\rm GSVD}^{(l)}\geq
0\end{array}\right.$;

11
$\boldsymbol{\sigma}^{(k,l+1)}=\boldsymbol{\sigma}^{(k,l)}+\alpha_{l}\mathbf{q}_{\rm
GSVD}^{(l)}$;

12

13 $\boldsymbol{\sigma}^{(k+1)}=\boldsymbol{\sigma}^{(k,\infty)}$

14

15 $\boldsymbol{\sigma}^{(k+1)}=\mbox{{vec}}(\mathbf{\Sigma}^{(k+1)})$;

16 $\boldsymbol{\xi}^{(k,1)}=\boldsymbol{\sigma}^{(k+1)}$;

17 for  _$l=1,2,\ldots$_ do

18 $\widetilde{\mathbf{u}}^{(l)}=D\boldsymbol{\xi}^{(k,l)}$;

19
$\mathbf{u}^{(l)}=\widetilde{\mathbf{u}}^{(l)}\left(1-\left(\frac{(\widetilde{\mathbf{u}}^{(j)})^{2}+\varepsilon^{2}}{\varepsilon^{2}}\right)^{q/2-1}\right)$
;

20
$\boldsymbol{\xi}^{(k,l+1)}={C}^{T}\left({I}+\eta\Lambda^{T}\Lambda\right)^{-1}{C}(\boldsymbol{\sigma}^{(k+1)}+\eta{D}^{T}\mathbf{u}^{(j)})$;

21 $\mbox{{vec}}(\mathbf{\Xi}^{(k+1)})=\boldsymbol{\xi}^{(k,\infty)}$;

22

23

Algorithm 1 Parallel alternating minimization for the solution of (9)

##  5 Numerical examples

In the numerical tests illustrated in this section we consider both synthetic
and experimental data. The two-dimensional representation of the electrical
conductivity is determined under the assumption that the magnetic permeability
is the one of the free space, i.e., $\boldsymbol{\mu}\equiv\mu_{0}$.

Throughout this section, we show the effectiveness of our method and we
compare the results obtained from the algorithm presented in this paper with
those obtained by applying the method described in [14], i.e., by solving each
one-dimensional problem independently and then “stacking” the obtained
solutions.

All the computations are performed on an Intel(R) Xeon(R) Gold 6136 CPU @
3.00GHz computer with 128Gb of RAM memory and 32 cores, running the Debian
GNU/Linux operating system and Matlab R2020b.

###  5.1 Synthetic data

In the first test we consider two different configurations of the FDEM device.
We show that the proposed algorithm is able to reconstruct an accurate
approximation of the electrical conductivity for both devices. In out second
test we report the solution obtained by varying the dimensions of the problem
discretization for one configuration.

We generate the data matrix $\mathbf{B}_{\rm{exact}}$ of a chosen dimension
$m\times N$ and we simulate the presence of noise in the data, by letting

$\mathbf{e}=\frac{\delta}{\sqrt{m}}\left\|\mathbf{B}_{\rm{exact}}\right\|_{F}^{2}\mathbf{w},$

where $\mathbf{w}$ is a vector with normally distributed entries having zero
mean and unit variance, and $\delta$ represents the noise level. The new data
matrix is denoted by $\mathbf{B}^{\delta}$.

#### Test 1.

This example concerns the reconstruction of a two-dimensional model generated
from $50$ soundings along a $10$ m straight line and $20$ layers in detph. It
is characterized by an increasing change in conductivity (from $0$ S/m to $1$
S/m) occurring at an increasing depth.

The synthetic data simulate an acquisition performed by both the Geophex GEM-2
and CMD Explorer instruments, with two orientations of the coils and one
measurement height $h=1$ m, i.e., the instrument is kept $1$ m above the
ground. We consider a slice of the ground that is $10$ meters wide. The first
device, the Geophex GEM-2, works with an intercoil distance of
$1.66\leavevmode\nobreak\ {\rm m}$ ($\rho=1.66\leavevmode\nobreak\ {\rm m}$)
and six different frequencies
$f=775,1175,3925,9825,21725,47025\leavevmode\nobreak\ {\rm Hz}$, while the
second one, the CDM Explorer, measures with three different intercoil
distances $\rho=1.48,2.82,4.49\leavevmode\nobreak\ {\rm m}$ and only one
frequency $f=10\leavevmode\nobreak\ {\rm kHz}$. The data values are finally
perturbed by uncorrelated Gaussian noise with standard deviation
$\delta=10^{-2}$. Here, the choice of the regularization parameter $\ell$ does
not affect substantially the computed results. This is due to the structure of
the Jacobian matrix defined in (18). The matrix $\widehat{J}$ is obtained by
stacking the Jacobian of ${M}$ and a scaled identity matrix, therefore, the
number of the singular values does not change, but they do not decrease as
rapidly as the one of $\widetilde{J}$. Therefore, we set in all our
experiments $\ell=15$. From our experience we observed that a larger value of
$\ell$ does not neither improve nor deteriorate the obtained approximate
solutions. In the algorithm described in [14] the Jacobian matrix of the
problem was very ill-conditioned and the results depended on the choice of the
truncation parameter. We have established a maximum number of $50$ iterations,
$p=2$, $q=0.1$, the regularization parameter $\gamma=10^{-4}$, and
$\sigma_{0}=0.1$.

In the first row of Figure 1, we report the exact solution used to generate
the data using both the devices: Geophex - GEM 2 on the left, CDM Explorer on
the right. On the middle row, we report the reconstruction of the one-
dimensional models side by side in a pseudo two-dimensional fashion as in
[14]. As we described above, the method in [14] solves independently the $N$
one-dimensional problems in (2). Then the solution is visualized by stacking
the one-dimensional reconstructions obtaining a two-dimensional image.
Finally, the obtained results by the algorithm presented in this paper are
depicted in the last row.

We can see that, for both the instruments configuration, the reconstruction of
the electrical conductivity obtained from the variational model described here
is much more accurate than the ones obtained by the method described in [14],
in which there is not lateral continuity in the results. This shows that the
additional regularization term in (3) significantly helps the reconstruction
of an approximate solution. In Figure 2 we report the reconstruction of a
single sounding (the 16th column of the two-dimensional synthetic model of
Figure 1 for both the devices). We compare the one-dimensional solution
computed by the algorithm in [14] to the one obtained with the proposed
method. We can observe that the introduction of the regularization term allows
us to improve the accuracy of the reconstruction. Finally, we compute the
Relative Restoration Error (RRE) for each computed solution, defined by

${\rm RRE}(\mathbf{\Sigma})=\frac{\left\|\mathbf{\Sigma}-\mathbf{\Sigma}_{\rm
exact}\right\|_{F}}{\left\|\mathbf{\Sigma}_{\rm exact}\right\|_{F}},$

where $\mathbf{\Sigma}_{\rm exact}$ denotes the exact solution of the problem.
We report the obtained results in Table 1. We can observe that the RRE
obtained with our proposal is significantly lower than the one obtained with
the method proposed in [14].

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1: Reconstruction of the electrical conductivity from two devices
configuration. The first row represents the exact solution used to generate
the data with both the Geophex - GEM 2 (panel (a)) and the CDM Explorer (panel
(b)). On the second row, is reported the reconstruction of the one-dimensional
models side by side in a pseudo two-dimensional fashion (panels (c) and (d)).
On the last row, the obtained results by the algorithm presented in this paper
are depicted (panels (e) and (f)). The data values are corrupted by Gaussian
noise with noise level $\delta=10^{-2}$, the maximum number of iterations is
fixed at $50$, $p=2$, $q=0.1$, and the regularization parameter
$\gamma=10^{-4}$.

(a)

(b)

Figure 2: One-dimensional reconstruction of the two-dimensional synthetic
model generated in Figure 1. On the left, $\boldsymbol{\sigma}_{16}$ in the
model generated by Geophex - GEM 2 (panel (a)). On the right
$\boldsymbol{\sigma}_{16}$ in the model generated by CDM Explorer (panel (b)).
The exact solution is depicted by a black dotted line, the profile computed by
the method in [14] is reported with a red dashed line, and the profile
obtained with our method is shown by a blue solid line.

#### Test 2.

This example reports the results for different dimensions of the matrix
solution: $50$, $100$ and $200$ soundings along a $10$ m straight line and
$20$, $50$ and $100$ layers, respectively. In this case, the synthetic data is
generated by the Geophex GEM-2 ($\rho=1.66\leavevmode\nobreak\ {\rm
m},\;f=775,\;1175,\;3925,\;9825,\;21725,\;47025\leavevmode\nobreak\ {\rm Hz}$)
with two orientations of the coils and one measurement height $h=1$ m. We use
the same parameters as in the previous example.

In the first column of Figure 3, we report the reconstructions of
$\mathbf{\Sigma}_{1}\in{\mathbb{R}}^{20\times 50}$. The second column depicts
the reconstructions of $\mathbf{\Sigma}_{2}\in{\mathbb{R}}^{50\times 100}$.
The last column increases the dimension of the exact solution to
$\mathbf{\Sigma}\in{\mathbb{R}}^{100\times 200}$. In all the experiments the
right-hand sides are affected by Gaussian noise with $\delta=10^{-2}$, the
maximum number of iterations is equal to $50$, $p=2$, $q=0.1$, the
regularization parameter $\gamma=10^{-4}$, and $\sigma_{0}=0.2$.

From the visual inspection of the reconstructions in Figure 3 we can observe
that the proposed model is able to capture more accurately the structure of
the soil in all proposed examples. Moreover, we can see that, since the
reconstruction is not obtained by simply stacking one-dimensional vectors, the
approximate solutions computed by our proposal are smoother than the ones
obtained with the method described in [14]. This is confirmed by the values of
the RRE reported in Table 1.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 3: Reconstruction of the electrical conductivity from data generated by the Geophex - GEM 2. The first row reports the exact solutions, the second row shows the reconstructions obtained with the method in [14], and the third row collects the approximate solution computed by our algorithmic proposal. Each column reports the results obtained with different examples. In the first column we consider $\mathbf{\Sigma}_{1}\in{\mathbb{R}}^{20\times 50}$, in the second column $\mathbf{\Sigma}_{2}\in{\mathbb{R}}^{50\times 100}$, and in the third column $\mathbf{\Sigma}_{3}\in{\mathbb{R}}^{100\times 200}$. The data values are affected by Gaussian noise with noise level $\delta=10^{-2}$, the maximum number of iterations is equal to $50$, $p=2$, $q=0.1$, and the regularization parameter $\gamma=10^{-4}$. Table 1: Relative Restoration Error (RRE) obtained in all the considered examples on synthetic data with the algorithm proposed in [14] and our proposal. Example | Method | RRE
---|---|---
Test 1 | Geophex | [14] | $0.55439$
Algorithm 1 | $0.37832$
Explorer | [14] | $0.41573$
Algorithm 1 | $0.35842$
Test 2 | $20\times 50$ | [14] | $0.28644$
Algorithm 1 | $0.25646$
$50\times 100$ | [14] | $0.65784$
Algorithm 1 | $0.36258$
$100\times 200$ | [14] | $0.93571$
Algorithm 1 | $0.36258$

###  5.2 Experimental data

In this section, we consider an experimental data set collected with a
multiconfiguration EMI device at the Molentargius Saline Regional Nature Park,
located east of Cagliari in southern Sardinia, Italy. The dataset was first
studied in [14]; see [14] for a description of the geographical location of
the test site.

As described in [14], the CDM Explorer
($\rho=1.48,2.82,4.49\leavevmode\nobreak\ {\rm m}$ and
$f=10\leavevmode\nobreak\ {\rm kHz}$) was used to collect the EMI data along a
$200$ m straight-line and was carried at $0.9$ m of height above the ground,
providing three simultaneous measurements for each orientation of the device.
The data were collected in continuous mode, with a $0.5$ s time step, first
using the horizontal orientation and then the vertical one. With the aim of
having the same number of equally spaced measurement points, the authors of
[14] merged the data by spatially resampling at $0.5$ m intervals from a
common starting point. This allowed to set up a dataset consisting of a series
of $401$ depth soundings with six complex (quadrature and in-phase components)
measurements each, suitable for recovering the soil electrical conductivities
in order to reconstuct the water table. The data set is available at the web
page http://bugs.unica.it/cana/datasets/.

(a)

(b)

Figure 4: Reconstruction of the electrical conductivity from real data
collected at Molentargius Saline Regional Nature Park in Cagliari by the CDM
Explorer. The panel on the left reports the reconstruction obtained with the
method in [14], and the right one collects the approximate solution computed
by our algorithmic proposal. In both panels we consider
$\mathbf{\Sigma}\in{\mathbb{R}}^{50\times 401}$ and a maximum number of
iterations equal to $10^{3}$. In panel (a) we use the first derivative
operator as the regularization matrix, while in panel (b) we fix $p=2$,
$q=0.1$, and we set the regularization parameter $\gamma=10^{-4}$.

In Figure 4, we compare the results obtained by the method described in [14]
where we used the first derivative operator as the regularization matrix, with
the ones computed by Algorithm 1 setting $p=2$, $q=0.1$, and the
regularization parameter $\gamma=10^{-4}$. In both cases the maximum number of
iterations has been fixed to $10^{3}$. We discretize the soil with $50$ layers
leading to a reconstruction $\mathbf{\Sigma}\in{\mathbb{R}}^{50\times 401}$.

From a graphic examination of the reconstructions in Figure 4, we can observe
that the proposed method (Figure 4(b)) is able to capture more accurately the
structure of the soil, avoiding any “splicing”. On the other hand the
reconstruction obtained with the algorithm in [14], shown in Figure 4(a), is
“spliced” in several points, that are highlighted in Figure 5. Moreover, we
can observe that both reconstructions have the same shape and structure,
showing the water table interface, even if it is not easy to identify its
exact depth.

We remark here that in these experiments, the same starting model was adopted
for both of the reconstructions but still, there is an important dependence on
the initial solution of the iterative method, also in our optimization
algorithm.

(a)

(b)

Figure 5: Plot of the first 100 soundings of the reconstructions in Figure 4.

##  6 Conclusions

In this paper we have proposed a non-linear model for the inversion of FDEM
data. Although the two-dimensional problem can be seen as a sequence of
“stacked” independent one-dimensional problems, the proposed model couples
them. This allows us to exploit the “horizontal” information to largely
improve the quality of the computed solution avoiding the “splicing” that
occurs when we deal with noisy data.

We were able to show that the outlined minimization problem has a solution,
albeit non-unique, and that it induces a regularization method. We provided an
algorithm for the solution of the problem and showed its performances on some
synthetic and real data. From these experiments we were able to show that our
method reliably computes accurate solutions even when compared with state-of-
the-art method like the one in [14].

Matters of future research include an efficient implementation exploiting
Krylov subspaces. Moreover, we plan on using the Alternating Direction
Multiplier Method (ADMM) and its accelerations (see, e.g., [4, 19, 6]) for
efficiently computing a solution of (3).

## Acknowledgements

The authors would like to thank Prof. G. Deidda for providing us with the
experimental data and, more in general, for his help, Prof. G. Rodriguez for
the precious discussions and the anonymous referees for their comments that
greatly improved the quality of this paper. The authors are members of the
GNCS group of INdAM. Alessandro Buccini is partially supported by Regione
Autonoma della Sardegna research project “Algorithms and Models for Imaging
Science [AMIS]” (RASSR57257, intervento finanziato con risorse FSC 2014-2020 -
Patto per lo Sviluppo della Regione Sardegna). Patricia Díaz de Alba is
partially supported by INdAM-GNCS 2020 project “Tecniche numeriche per
l’analisi delle reti complesse e lo studio dei problemi inversi”.

## References

  * [1] Z.-Z. Bai, A. Buccini, K. Hayami, L. Reichel, J.-F. Yin, and N. Zheng, Modulus-based iterative methods for constrained Tikhonov regularization, Journal of Computational and Applied Mathematics, 319 (2017), pp. 1–13.
  * [2] Å. Björck, Numerical Methods for Least Squares Problems, SIAM, Philadelphia, 1996.
  * [3] B. Borchers, T. Uram, and J. Hendrickx, Tikhonov regularization of electrical conductivity depth profiles in field soils, Soil Science Society of America Journal, 61 (1997), pp. 1004–1009.
  * [4] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends in Machine Learning, 3 (2011), pp. 1–122.
  * [5] A. Buccini, O. De la Cruz Cabrera, M. Donatelli, A. Martinelli, and L. Reichel, Large-scale regression with non-convex loss and penalty, Applied Numerical Mathematics, 157 (2020), pp. 590–601.
  * [6] A. Buccini, P. Dell’Acqua, and M. Donatelli, A general framework for ADMM acceleration, Numerical Algorithms, 85 (2020), pp. 829–848.
  * [7] A. Buccini, M. Pasha, and L. Reichel, Modulus-based iterative methods for constrained $\ell_{p}$-$\ell_{q}$ minimization, Inverse Problems, 36 (2020), p. 084001.
  * [8] A. Buccini and L. Reichel, An $\ell^{2}-\ell^{q}$ regularization method for large discrete ill-posed problems, Journal of Scientific Computing, 78 (2019), pp. 1526–1549.
  * [9] , An $\ell^{p}-\ell^{q}$ minimization method with cross-validation for the restoration of impulse noise contaminated images, Journal of Computational and Applied Mathematics, (2020), p. 112824.
  * [10] , Generalized cross validation for $\ell^{p}-\ell^{q}$ minimization, Numerical Algorithms, In press (2021), pp. 1–22.
  * [11] G. P. Deidda, P. Díaz de Alba, C. Fenu, G. Lovicu, and G. Rodriguez, FDEMtools: a MATLAB package for FDEM data inversion, Numerical Algorithms, 84(4) (2020), pp. 1313–1327.
  * [12] G. P. Deidda, P. Díaz de Alba, and G. Rodriguez, Identifying the magnetic permeability in multi-frequency EM data inversion, Electronic Transactions on Numerical Analysis, 47 (2017), pp. 1–17.
  * [13] G. P. Deidda, P. Díaz de Alba, G. Rodriguez, and G. Vignoli, Smooth and sparse inversion of EMI data from multi-configuration measurements, in 2018 IEEE 4th International Forum on Research and Technology for Society and Industry (RTSI) (RTSI 2018), Palermo, Italy, Sept. 2018, pp. 213–218.
  * [14] , Inversion of Multiconfiguration Complex EMI Data with Minimum Gradient Support Regularization: A Case Study, Mathematical Geosciences, 52(7) (2020), pp. 945–970.
  * [15] G. P. Deidda, C. Fenu, and G. Rodriguez, Regularized solution of a nonlinear problem in electromagnetic sounding, Inverse Problems, 30 (2014), p. 125014 (27 pages).
  * [16] P. Díaz de Alba, L. Fermo, C. van der Mee, and G. Rodriguez, Recovering the electrical conductivity of the soil via a linear integral model, Journal of Computational and Applied Mathematics, 352 (2019), pp. 132–145.
  * [17] P. Díaz de Alba and G. Rodriguez, Regularized inversion of multi-frequency EM data in geophysical applications, in Trends in Differential Equations and Applications, F. Ortegón Gallego, M. Redondo Neble, and J. Rodríguez Galván, eds., vol. 8 of SEMA SIMAI Springer Series, Springer, Switzerland, 2016, pp. 357–369.
  * [18] H. Engl, M. Hanke, and A. Neubauer, Regularization of Inverse Problems, Kluwer, Dordrecht, 1996.
  * [19] T. Goldstein, B. O’Donoghue, S. Setzer, and R. Baraniuk, Fast alternating direction optimization methods, SIAM Journal on Imaging Sciences, 7 (2014), pp. 1588–1623.
  * [20] G. H. Golub and C. F. Van Loan, Matrix Computations, The John Hopkins University Press, Baltimore, third ed., 1996.
  * [21] L. Grippof and M. Sciandrone, Globally convergent block-coordinate techniques for unconstrained optimization, Optimization methods and software, 10 (1999), pp. 587–637.
  * [22] P. C. Hansen, Rank–Deficient and Discrete Ill–Posed Problems, SIAM, Philadelphia, 1998.
  * [23] J. M. H. Hendrickx, B. Borchers, D. L. Corwin, S. M. Lesch, A. C. Hilgendorf, and J. Schlue, Inversion of soil conductivity profiles from electromagnetic induction measurements, Soil Science Society of America Journal, 66 (2002), pp. 673–685.  Package NONLINEM38 available at http://infohost.nmt.edu/~borchers/nonlinem38.html.
  * [24] B. Hofmann, B. Kaltenbacher, C. Poeschl, and O. Scherzer, A convergence rates result for Tikhonov regularization in Banach spaces with non-smooth operators, Inverse Problems, 23 (2007), pp. 987–1010.
  * [25] G. Huang, A. Lanza, S. Morigi, L. Reichel, and F. Sgallari, Majorization–minimization generalized Krylov subspace methods for $\ell_{p}-\ell_{q}$ optimization applied to image restoration, BIT Numerical Mathematics, 57 (2017), pp. 351–378.
  * [26] J. D. McNeill, Electromagnetic terrain conductivity measurement at low induction numbers, Technical Report TN-6 Geonics Limited, 1980.
  * [27] M. K. Ng, R. H. Chan, and W.-C. Tang, A fast algorithm for deblurring models with neumann boundary conditions, SIAM Journal on Scientific Computing, 21 (1999), pp. 851–866.
  * [28] L. Reichel and G. Rodriguez, Old and new parameter choice rules for discrete ill-posed problems, Numerical Algorithms, 63 (2013), pp. 65–87.
  * [29] L. I. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D: nonlinear phenomena, 60 (1992), pp. 259–268.
  * [30] A. Tikhonov and D. Shakhsuvarov, On the possibility of using the impedance of the natural electromagnetic field of the earth in exploration of the upper layers, Izvest. Akad. Nauk SSSR, Set. Geofiz, (1956), pp. 410–418.
  * [31] J. R. Wait, Geo–Electromagnetism, New York: Academic, New York, NY, 1982.