# Sparse recovery for inverse potential problems in divergence form

Laurent Baratchart, Cristóbal Villalobos Guillén, Douglas Hardin, Juliette Leblond

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

## To cite this version:

Laurent Baratchart, Cristóbal Villalobos Guillén, Douglas Hardin, Juliette Leblond. Sparse recovery for inverse potential problems in divergence form. 9th International Conference on New Computational Methods for Inverse Problems, NCMIP 2019 24 May 2019, Cachan, France, 1476, IOP Publishing, pp.012009, 2020, Journal of Physics: Conference Series, 10.1088/1742-6596/1476/1/012009. hal-03434756

HAL Id: hal-03434756

https://inria.hal.science/hal-03434756

Submitted on 18 Nov 2021

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

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

<figure><img src="image_2.png" /></figure>
# Sparse recovery for inverse potential problems in divergence form

Laurent Baratchart, Cristobal Villalobos-Guillen, Douglas Hardin, Juliette Leblond

INRIA, 2004 route des lucioles, 06902 Sophia Antipolis cedex, France

E-mail: Laurent.Baratchart@inria.fr

**Abstract.** We discuss recent results from [10] on sparse recovery for inverse potential problem with source term in divergence form. The notion of sparsity which is set forth is measure-theoretic, namely pure 1-unrectifiability of the support. The theory applies when a superset of the support is known to be slender, meaning it has measure zero and all connected components of its complement has infinite measure in $\mathbb{R}^3$. We also discuss open issues in the non-slender case.

## 1. Introduction

Inverse potential problems with source term in divergence form consist in recovering a $\mathbb{R}^3$-valued distribution $\mu$, knowing the potential $\Phi$ of $\operatorname{div}\mu$ which is the solution to the Poisson-Hodge equation $\Delta\Phi = \operatorname{div}\mu$ having “least growth” at infinity. In practice, a superset $S$ of the support of $\mu$ is known a priori, and sensors will measure the field $\nabla\Phi$ rather than the potential $\Phi$ itself.

Issues of this kind typically arise in source identification from field measurements for Maxwell's equations, in the quasi-static regime. They occur for instance in electroencephalography (EEG), magneto-encephalography (MEG), geomagnetism and paleomagnetism, as well as in several non-destructive testing problems, see e.g. [1, 2, 3, 4, 5] and their bibliographies. A model problem of our particular interest is inverse scanning magnetic microscopy, as considered for instance in [9, 6, 7, 8] to recover magnetization distributions of thin rock samples, but the considerations below are of a rather general and mathematical nature.

Such problems are known to be difficult, for they are not only ill-posed but the forward operator, mapping $\mu$ to the field, is not even injective. Recently, in the preprint [10], notions of sparsity have been introduced concerning $\mu$, when the latter is a finite $\mathbb{R}^3$-valued measure. They justify the use of Tikhonov-like regularization schemes that minimize the residuals while penalizing the total variation norm, in order to asymptotically reconstruct a sparse measure $\mu$ when the regularization parameter goes to zero, under a specific assumption on $S$: it should be *slender*, meaning it has measure zero and each connected component of $\mathbb{R}^3 \setminus S$ has infinite measure.

This situation is reminiscent of compressive sensing, where sparse solutions to underdetermined systems of linear equations in $\mathbb{R}^n$ (i.e. solutions having a large number of zero components) are sought by minimizing the residuals while penalizing the $l^1$-norm; the gist of this approach is that, for “most” large underdetermined systems, the solution with minimal $l^1$-norm is also the sparsest solution, see e.g. [11].
However, in the present, infinite-dimensional context, it is unclear which assumption on $\mu$ will ensure that it has minimum total variation among all solutions to the (noise-free) inverse problem, and why such an assumption should connect with some kind of sparsity. In fact, the answer to such questions will much depend on the null-space of the forward operator. In [10], the assumption that $S$ is slender is to the effect that this kernel consists exactly of divergence-free $\mathbb{R}^2$-valued measures, also known as solenoids. Then, using a characterization of solenoids as integrals of elementary ones supported on curves [12], a natural notion of sparsity is found that ensures a sparse measure is mutually singular to all solenoids. This notion of sparsity which pertains to geometric measure theory, namely the support of $\mu$ should be purely 1-unrectifiable; roughly speaking this means it contains no rectifiable arc. For instance, a countable sum of Dirac masses will satisfy this condition, but other, more complicated supports would also qualify.

The goal of this paper is to present main results from [10], and to discuss new issues arising when $S$ is not slender.

We mention that a general Tikhonov-like regularization theory was developed in [14, 15, 16] for linear equations whose unknown is a $\mathbb{R}^n$-valued measure, by minimizing the residuals while penalizing the total variation. As expected from the non-reflexive character of spaces of measures, consistency estimates hold in the sense of weak-* convergence of subsequences to solutions of minimum total variation, or convergence in the Bregman distance when the so-called source condition holds (which is, by the way, not the case here). In principle, such methods yield algorithms to approximate a solution of minimum total variation to the initial equation by a sequence of discrete measures, but imply nothing about the nature of limit points as discrete measures are weak-* dense in the space of measures supported on an open subset of $\mathbb{R}^n$.

We note also that an infinite-dimensional recovery result for sparse measures, in the sense of being a sum of Dirac masses, was established in [17] for 1-D deconvolution issues, where a train of spikes is to be recovered from filtered observation thereof. Thus, [10] does not state the first sparse recovery result in an infinite-dimensional setting. It seems hoewever, that [10] produces the first sparse recovery result for convolution operators in space-dimension greater than 1. Moreover, we should stress in our case that the convolution kernel is singular and the null-space of the forward operator has infinite dimension.

## 2. The inverse problem

Without loss of generality, we consider the issue of recovering a magnetization distribution from a collection of measurements of the magnetic field the magnetization generates. For a closed subset $S \subset \mathbb{R}^3$, let $\mathcal{M}(S)$ denote the space of finite signed Borel measures on $\mathbb{R}^3$ whose support lies in $S$. We model magnetization distributions supported in $S$ as $\mathbb{R}^3$-valued measures $\mu \in \mathcal{M}(S)^3$. Hereafter, we often call a member of $\mathcal{M}(S)^3$ a magnetization supported on $S$, as this terminology is suggestive of the problems we adress.

The magnetic field $b(\mu)$ generated by a magnetization $\mu \in \mathcal{M}(S)^3$ is defined, at a point $x$ not in the support of $\mu$, by the formula [18]:

$$b(\mu)(x) := -c \left( \int \frac{1}{|x-y|^3} d\mu(y) - 3 \int (x-y) \frac{(x-y) \cdot d\mu(y)}{|x-y|^5} \right) \quad x \notin \text{supp } \mu, \quad (1)$$

where $c = 10^{-7}Hm^{-1}$ and for $x, y \in \mathbb{R}^3$, $x \cdot y$ and $|x|$ denote the Euclidean scalar product and norm. Equivalently:

$$b(\mu)(x) = -\mu_0 \nabla \Phi(\mu)(x), \quad x \notin \text{supp } \mu, \quad (2)$$

with $\mu_0 = 4\pi \times c$ and $\Phi(\mu)$ is the scalar magnetic potential defined by

$$\Phi(\mu)(x) := \frac{1}{4\pi} \int \nabla_y \frac{1}{|x-y|} \cdot d\mu(y) = \frac{1}{4\pi} \int \frac{x-y}{|x-y|^3} \cdot d\mu(y). \quad (3)$$
and $\nabla_y$ denotes the gradient with respect to $y$

Generally speaking, the inverse magnetization problem is to recover $\mu \in \mathcal{M}(S)^3$, knowing $b(\mu)$ in $\mathbb{R}^3 \setminus S$. However $b(\mu)$ is usually measured on a rather small subset $Q \subset \mathbb{R}^3 \setminus S$, typically a compact surface patch. Also, in most cases, only one component of $b(\mu)$ can be measured, because coils are oriented. For the sake of simplicity, we shall assume that $S$ has connected complement and is contained in the closed lower half-space $H := \{x = (x_1, x_2, x_3)^T \in \mathbb{R}^3, x_3 \le 0\}$, while $Q$ is a compact set of Hausdorff dimension greater than 1, contained in a horizontal plane $\Pi_h$ at strictly positive height $h$. Also, the component of $b(\mu)$ which is measured will be $b_3(\mu)$, the third (vertical) one. This is the setting adopted in scanning magnetic microscopy [9, 6, 7, 8]. One could also consider the case where $S$ is a surface disconnecting the space (e.g. a plane), in which case $Q$ should intersect each of component of $\mathbb{R}^3 \setminus S$ and be contained in a union of real analytic surfaces positively separated from $S$ and satisfying mild conditions. We refer to [10] for this more exhaustive framework.

Letting $\{e_j, 1 \le j \le 3\}$ indicate the canonical basis of $\mathbb{R}^3$, we get from (1) that

$$b_3(\mu)(x) := -\frac{\mu_0}{4\pi} \int \mathbf{K}(x-y) \cdot d\mu(y), \qquad (4)$$

where

$$K(x) = \frac{e_3}{|x|^3} - 3x \frac{x_3}{|x|^5} = \nabla \left( \frac{x_3}{|x|^3} \right). \qquad (5)$$

We define the forward operator $A: \mathcal{M}(S)^3 \to L^2(Q)$ by

$$A(\mu)(x) := b_3(\mu)(x), \quad x \in Q. \qquad (6)$$

Now, the inverse problem consists in recovering $\mu$ knowing $A(\mu)$.

Note that in practice, only pointwise values of $b_3(\mu)$ can be estimated, whereas we assume here knowledge of $b_3(\mu)$ at each point of $Q$. We shall ignore this important issue, as it pertains to a numerical approach of the the problem which is beyond the scope of the present paper, devoted to basic principles.

## 2.1. Slenderness and null-space of the forward operator

Let $\mathcal{L}_3$ denote Lebesgue measure on $\mathbb{R}^3$. We say that $E \subset \mathbb{R}^3$ is slender if $\mathcal{L}_3(E) = 0$ and any connected component $C$ of $\mathbb{R}^3 \setminus E$ has $\mathcal{L}_3(C) = +\infty$.

Clearly, the potential $\Phi(\mu)$ defined by (3) is harmonic in $\mathbb{R}^3 \setminus S$, and so are the components of $b(\mu)$. It is easy to check that $\Phi(\mu)$ extends to a function in $L^2(\mathbb{R}^3) + L^1(\mathbb{R}^3)$ (still denoted $\Phi(\mu)$), and that $b(\mu)$ extends to a $\mathbb{R}^3$-valued divergence-free distribution [10, Prop. 2.1], with

$$\Delta\Phi = \operatorname{div}\mu \quad \text{and} \quad b(\mu) = \mu_0(\mu - \nabla\Phi(\mu)). \qquad (7)$$

It is not difficult to check that $A(\mu)$ characterizes $b(\mu)$ completely, and we explain this in the simple case where $S$ has connected complement: if $A(\mu) = 0$, then $b_3(\mu) = 0$ on $Q$ and consequently on the entire plane $\Pi_h$, because it is real analytic in $\{x_3 > 0\}$, being harmonic in the upper half-space. Then, $\Phi$ is a harmonic function in the upper half-space which solves a Neumann problem in $\{x_3 > h\}$ with vanishing normal derivative, hence it is constant and since it lies in $L^2(\mathbb{R}^3) + L^1(\mathbb{R}^3)$ it must be zero.

Still, $A$ will generally have nontrivial null-space, because the mapping $\mu \to b(\mu)$ is typically not injective. We say that $\mu, \nu \in \mathcal{M}(S)^3$ are S-equivalent if $b(\mu)$ and $b(\nu)$ agree on $\mathbb{R}^3 \setminus S$, in which case we write $\mu \equiv \nu[S]$. A magnetization $\mu$ is said to be S-silent if $\mu \equiv 0[S]$; i.e., if $b(\mu)$ vanishes on $\mathbb{R}^3 \setminus S$.
Since no nonzero harmonic function can lie in $L^2(\mathbb{R}^3) + L^1(\mathbb{R}^3)$, by the mean value property and Liouville's theorem, it follows from (7) that a divergence $\mu$ belongs to the kernel of $A$. The converse needs not hold in general, but it does if $S$ is slender:

**Theorem.** If $S$ is a slender set and $\mu$ is $S$-silent, then $\mathrm{div}\,\mu = 0$.

For a proof, we refer to [10, Thm. 2.2].

**Corollary.** If $S$ is slender, the null-space of $A$ consists of all divergence-free $\mathbb{R}^3$-valued measures on $\mathbb{R}^3$ that are supported on $S$.

## 2.2. Divergence-free measures, pure 1-unrectifiability and total variation

We let $\mathcal{H}_1$ indicate 1-dimensional Hausdorff measure, see [13] for a definition. A set $E \subset \mathbb{R}^3$ is said to be 1-rectifiable if there exist Lipschitz maps $f_i: \mathbb{R} \to \mathbb{R}^3$, $i = 1, 2, \dots$, such that

$$ \mathcal{H}_1 \left( E \setminus \bigcup_{i=1}^{\infty} f_i(\mathbb{R}) \right) = 0. $$

A set $B \subset \mathbb{R}^3$ is purely 1-unrectifiable if $\mathcal{H}_1(E \cap B) = 0$ for every 1-rectifiable set $E$, see [19] for these definitions. Clearly a set of $\mathcal{H}_1$-measure zero is purely 1-unrectifiable.

For $\mu \in \mathcal{M}(\mathbb{R}^3)$ the total variation measure $|\mu|$ is defined on Borel sets $B \subset \mathbb{R}^3$ by

$$ |\mu|(B) := \sup_{\mathcal{P}} \sum_{P \in \mathcal{P}} |\mu(P)|, \qquad (8) $$

where the supremum is taken over all finite Borel partitions $\mathcal{P}$ of $B$. The total variation norm $\|\mu\|_{TV}$ is then defined as $|\mu|(\mathbb{R}^3)$.

The theorem below, which is crucial to the present approach of the inverse problem, entails that a magnetization with purely 1-unrectifiable support is the unique element of minimal total variation norm in its coset modulo the null-space of $A$.

**Theorem** Suppose $S$ is slender. If $\mu \in \mathcal{M}(S)^3$ has support that is purely 1-unrectifiable and $\nu \in \mathcal{M}(S)^3$ is $S$-equivalent to $\mu$, then $\|\nu\|_{TV} > \|\mu\|_{TV}$ unless $\nu = \mu$.

The proof rests on [12, Thm. A] which represents divergence-free measures as integrals of elementary measures of the form $R_\gamma(B) = \int_B \tau d(\mathcal{H}_1\lfloor\gamma)$, where $\gamma$ is an oriented Lipschitz arc and $\tau$ its unit tangent, with $\mathcal{H}_1\lfloor\gamma$ to mean the restriction of $\mathcal{H}_1$ to the image of $\gamma$.

## 3. Total variation regularization and consistency of sparse recovery

For $\mu \in \mathcal{M}(S)^3$, $f \in L^2(Q)$, and $\lambda > 0$, define

$$ \mathcal{F}_{f,\lambda}(\mu) := \|f - A\mu\|_{L^2(Q)}^2 + \lambda \|\mu\|_{TV}. \qquad (9) $$

To recover a magnetization $\mu_0$ from measurements $f$ of $A\mu_0$, we pick $\lambda > 0$ and consider the following regularization scheme: to find $\mu_\lambda \in \mathcal{M}(S)^3$ such that

$$ \mathcal{F}_{f,\lambda}(\mu_\lambda) = \inf_{\mu \in \mathcal{M}(S)^3} \mathcal{F}_{f,\lambda}(\mu). \qquad (10) $$

More precisely, to account for measurement noise, we assume that $f = f_e = A\mu_0 + e$ and we call $\mu_{\lambda,e}$ a minimizer of (10). It is easy to see that such a minimizer exists, and it follows from [14, Thms. 2&5] or [15, Thm. 3.5&4.4] that any weak-* limit point of the family $\mu_{\lambda,e}$ as both $\lambda$ and $\|e\lambda^{-1/2}\|_{L^2(Q)}$ tend to 0 is a magnetization $\nu$ such that $A\nu = A\mu_0$ of minimum total variation under this condition. In particular, if there is a unique such magnetization, we can easily formulate a weak-* sequential consistency result in the zero-noise limit, by letting $\lambda_n$ go to zero more slowly than $e_n\|_{L^2(Q)}^2$ (the so-called Morozov discrepancy principle).
The theorem below dwells on this and on the previous theorem, but goes a little further in that the convergence holds not only for $\mu_{\lambda,e}$ but also for $|\mu_{\lambda,e}|$, which is important for recovery (think of an oscillating density like $e^i n\theta$ on the unit circle, which goes weak-* to 0 as $n \to \infty$ but still has total variation $2\pi$ for each $n$). Moreover, convergence takes place in the narrow sense (test functions should be continuous and bounded but need not have compact support)

**Theorem** Let $S$ be slender and $\mu_0 \in \mathcal{M}(S)^3$ have purely 1-unrectifiable support. Then, $\mu_{\lambda,e}$ converges narrowly sequentially to $\mu_0$ and $|\mu_{\lambda,e}|$ converges narrowly sequentially to $|\mu_0|$ as $\lambda \to 0$ and $\|e\|_{L^2(Q)}/\sqrt{\lambda} \to 0$.

### 4. Issues in the non slender case

A typical slender set is two-dimensional. It could be a piece of plane, or the entire plane, or it could be a piece of sphere but not the entire sphere, nor a ball.

If $S$ is a genuine 3-D object, like the interior of a compact surface $\Sigma$, then it is not slender and the previous analysis fails. In fact, there are in this case magnetizations which are $S$-silent but not divergence-free: an example when $\Sigma$ is Lipschitz is given by the gradient of a function of bounded variation inside $\Sigma$ which has constant trace on $\Sigma$. Still, such a magnetization turns out to be singular with measures with pure 1-unrectifiable support. It would be most interesting to describe all silent magnetizations supported inside $\Sigma$. In this connection, we mention that such magnetizations, if they have $L^p$ density, must be the sum of a divergence-free field and a gradient as above. For finite measures, however, no characterization is known.

[1] S. Baillet, J.C. Mosher and R. M. Leahy 2001, Electromagnetic brain mapping, *IEEE Signal Processing Magazine*.

[2] H. T. Banks and F. Kojima 2002, Identification of material damage in two-dimensional domains using the SQUID-based nondestructive evaluation system, *Inverse Problems* **18** 1831-55.

[3] R. J. Blakely 1995, *Potential Theory in Gravity and Magnetic Applications*, Cambridge University Press.

[4] R. L. Parker 1994, *Geophysical inverse theory*, Princeton University Press.

[5] R. Kress, L. Kühn and R. Potthast 2002, Reconstruction of a current distribution from its magnetic field, *Inverse Problems* **18**, 1127-46.

[6] L. Baratchart and D.P. Hardin and E.A. Lima and E.B. Saff and B.P. Weiss 2013, Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions, *Inverse Problems* **29**, 15004.

[7] E. A. Lima, B. P. Weiss, L. Baratchart, D. P. Hardin and E. B. Saff 2013, Fast inversion of magnetic field maps of unidirectional planar geological magnetization, *Journal of Geophysical Research: Solid Earth* **118**, 2723-52.

[8] L. Baratchart and S. Chevillard and J. Leblond 2017, Silent and equivalent magnetic distributions on thin plates, *Harmonic Analysis, Function Theory, Operator Theory, and Their Applications, Theta series in advanced mathematics* **18**.

[9] B. P. Weiss, E. A. Lima, L. E. Fong, F. J. Baudenbacher 2007, Paleomagnetic analysis using SQUID microscopy, *Journal of Geophysical Research: Solid Earth* **112**.

[10] L. Baratchart, C. Villalobos-Guillen, D. P. Hardin, M. C. Northington and E.B. Saff 2018, Inverse Potential Problems for Divergence of Measures with Total Variation Regularization, *Preprint* https://arxiv.org/abs/1809.08334.

[11] S. Foucart and H. Rauhut 2013, *A Mathematical Introduction to Compressive Sensing*, Birkhäuser.

[12] S. K. Smirnov 1994, Decomposition of solenoidal vector charges into elementary solenoids, and the structure of normal one-dimensional flows, *St. Petersburg Math. J.* **5**, 841-67.

[13] L. C. Evans and R. F. Gariepy 1992, *Measure theory and fine properties of functions*, Studies in Advanced Mathematics, CRC Press.

[14] M. Burger and S. Osher 2004, Convergence rates of convex variational regularization, *Inverse Problems* **20**, 1411-21.

[15] B. Hoffmann, B. Kaltenbacher, C. Pöschl and O. Scherzer 2007, *Inverse Problems* **23**, 987-10.

[16] K. Bredies and H. K. Pikkarainen 2013, Inverse problems in spaces of measures, *ESAIM COCV* **19**, 190-08.

[17] V. Duval and G. Peyré 2015, Exact Support Recovery for Sparse Spikes Deconvolution, *FoCM* **15**, 1315-55.

[18] J. D. Jackson 1975, Classical electrodynamics, John Wiley & Sons.

[19] P. Mattila 1995, *Geometry of sets and measures in Euclidean spaces*, Cambridge University Press.