Imaging with highly incomplete and corrupted data

TL;DR

This paper investigates the limitations of $l_1$-minimization for imaging sparse scenes with noisy, incomplete data and proposes an augmented system with a noise collector to enhance robustness.

Contribution

It introduces a novel augmented linear system approach with a noise collector to improve $l_1$-minimization robustness in noisy imaging scenarios.

Findings

Standard $l_1$-minimization is not robust under high noise and resolution.

Augmented system with noise collector improves noise robustness.

Numerical results show effectiveness with noise collector size around 10 times the unknowns.

Abstract

We consider the problem of imaging sparse scenes from a few noisy data using an $l_1$-minimization approach. This problem can be cast as a linear system of the form $A \, \rho =b$, where $A$ is an $N\times K$ measurement matrix. We assume that the dimension of the unknown sparse vector $\rho \in {\mathbb{C}}^K$ is much larger than the dimension of the data vector $b \in {\mathbb{C}}^N$, i.e, $K \gg N$. We provide a theoretical framework that allows us to examine under what conditions the $\ell_1$-minimization problem admits a solution that is close to the exact one in the presence of noise. Our analysis shows that $l_1$-minimization is not robust for imaging with noisy data when high resolution is required. To improve the performance of $l_1$-minimization we propose to solve instead the augmented linear system $ [A \, | \, C] \rho =b$, where the $N \times \Sigma$ matrix $C$ is a noise collector. It is constructed so as its column vectors provide a frame on which the noise of the data, a vector of dimension $N$, can be well approximated. Theoretically, the dimension $\Sigma$ of the noise collector should be $e^N$ which would make its use not practical. However, our numerical results illustrate that robust results in the presence of noise can be obtained with a large enough number of columns $\Sigma \approx 10 K$.

Equations

$${\cal A}\,\mbox{\boldmath$\rho$}=\mbox{\boldmath$b$}\,,$$

$$\|\mbox{\boldmath{$\rho$}}\|_{\ell_{1}}\leqslant\gamma\,\|\mbox{\boldmath{$b$}}\|_{\ell_{2}}.$$

$$\mbox{\boldmath{$\rho$}}=[\rho_{1},\ldots,\rho_{K}]^{\intercal}\in\mathbb{C}^{K}\,,$$

$$\rho_{k}=\left\{\begin{array}[]{ll}\alpha_{j},&\hbox{ if }\|{\vec{\mbox{\boldmath{$z$}}}}_{j}-\vec{\mbox{\boldmath{$y$}}}_{k}\|_{\infty}<\hbox{ grid-size, for some }j=1,\ldots,M,\\ 0,&\hbox{ otherwise.}\end{array}\right.$$

$$\mbox{\boldmath{$g$}}(\vec{\mbox{\boldmath{$y$}}};\omega)=[G(\vec{\mbox{\boldmath{$x$}}}_{1},\vec{\mbox{\boldmath{$y$}}};\omega),G(\vec{\mbox{\boldmath{$x$}}}_{2},\vec{\mbox{\boldmath{$y$}}};\omega),\ldots,G(\vec{\mbox{\boldmath{$x$}}}_{N},\vec{\mbox{\boldmath{$y$}}};\omega)]^{\intercal}\,$$

$$G(\vec{\mbox{\boldmath{$x$}}},\vec{\mbox{\boldmath{$y$}}};\omega)=\frac{\exp\left(i\frac{\omega|\vec{\mbox{\boldmath{$x$}}}-\vec{\mbox{\boldmath{$y$}}}|}{c_{0}}\right)}{4\pi|\vec{\mbox{\boldmath{$x$}}}-\vec{\mbox{\boldmath{$y$}}}|}.$$

$$b(\vec{\mbox{\boldmath{$x$}}}_{r},\omega_{l})=\sum_{j=1}^{M}\alpha_{j}G(\vec{\mbox{\boldmath{$x$}}}_{r},{\vec{\mbox{\boldmath{$z$}}}}_{j};\omega_{l}).$$

$$\mbox{\boldmath{$b$}}=\frac{1}{\sqrt{NS}}[b(\vec{\mbox{\boldmath{$x$}}}_{1},\omega_{1}),b(\vec{\mbox{\boldmath{$x$}}}_{2},\omega_{1}),\dots,b(\vec{\mbox{\boldmath{$x$}}}_{N},\omega_{S})]^{\intercal}\,,$$

$${\cal A}=\frac{1}{\sqrt{NS}}\left(\begin{array}[]{cccc}\uparrow&\uparrow&&\uparrow\\ \mbox{\boldmath{$g$}}(\vec{\mbox{\boldmath{$y$}}}_{1};\omega_{1})&\mbox{\boldmath{$g$}}(\vec{\mbox{\boldmath{$y$}}}_{2};\omega_{1})&\ldots&\mbox{\boldmath{$g$}}(\vec{\mbox{\boldmath{$y$}}}_{K};\omega_{1})\\ \downarrow&\downarrow&&\downarrow\\ \uparrow&\uparrow&&\uparrow\\ \mbox{\boldmath{$g$}}(\vec{\mbox{\boldmath{$y$}}}_{1};\omega_{2})&\mbox{\boldmath{$g$}}(\vec{\mbox{\boldmath{$y$}}}_{2};\omega_{2})&\ldots&\mbox{\boldmath{$g$}}(\vec{\mbox{\boldmath{$y$}}}_{K};\omega_{2})\\ \downarrow&\downarrow&&\downarrow\\ \vdots&\vdots&&\vdots\\ \uparrow&\uparrow&&\uparrow\\ \mbox{\boldmath{$g$}}(\vec{\mbox{\boldmath{$y$}}}_{1};\omega_{S})&\mbox{\boldmath{$g$}}(\vec{\mbox{\boldmath{$y$}}}_{2};\omega_{S})&\ldots&\mbox{\boldmath{$g$}}(\vec{\mbox{\boldmath{$y$}}}_{K};\omega_{S})\\ \downarrow&\downarrow&&\downarrow\end{array}\right):=\left(\begin{array}[]{cccc}\uparrow&\uparrow&&\uparrow\\ \mbox{\boldmath{$a$}}_{1}&\mbox{\boldmath{$a$}}_{2}&..&\mbox{\boldmath{$a$}}_{K}\\ \downarrow&\downarrow&&\downarrow\end{array}\right).$$

$$\mbox{\boldmath$\rho$}_{\ell_{1}}=\mathop{\mbox{argmin}}\|\mbox{\boldmath{$\rho$}}\|_{\ell_{1}},\hbox{ subject to }{\cal A}\mbox{\boldmath{$\rho$}}=\mbox{\boldmath{$b$}}.$$

$$F(\mbox{\boldmath{$\rho$}},\mbox{\boldmath{$z$}})=\tau\|\mbox{\boldmath{$\rho$}}\|_{\ell_{1}}+\frac{1}{2}\|{\mathcal{A}}\mbox{\boldmath{$\rho$}}-\mbox{\boldmath$b$}\|^{2}_{\ell_{2}}+\langle\mbox{\boldmath{$z$}},\mbox{\boldmath$b$}-{\mathcal{A}}\mbox{\boldmath{$\rho$}}\rangle$$

$$\max_{\mbox{\boldmath{$z$}}}\min_{\mbox{\boldmath{$\rho$}}}F(\mbox{\boldmath{$\rho$}},\mbox{\boldmath{$z$}}).$$

$$|\langle\mbox{\boldmath{$a$}}_{i},\mbox{\boldmath{$a$}}_{j}\rangle|<\frac{1}{2M},\quad\forall i\neq j,$$

$$S_{j}=\left\{k~{}:~{}|\langle\mbox{\boldmath{$a$}}_{k},\mbox{\boldmath{$a$}}_{j}\rangle|\geqslant\frac{1}{3M}\right\}.$$

$${\bf Co}(\mbox{\boldmath{$\rho$}},\mbox{\boldmath{$\eta$}})=\sum_{j\in T}\left|\rho_{j}-\sum_{k\in S_{j}}\langle\mbox{\boldmath{$a$}}_{j},\mbox{\boldmath{$a$}}_{k}\rangle\eta_{k}\right|,$$

$${\bf In}(\mbox{\boldmath{$\rho$}},\mbox{\boldmath{$\eta$}})=\sum_{k\not\in\Upsilon}|\eta_{k}|,~{}\,\,~{}\Upsilon=\cup_{j\in T}S_{j}.$$

$$\gamma=\sup_{\mbox{\boldmath{$c$}}}\frac{\|\mbox{\boldmath{$\xi$}}\|_{\ell_{1}}}{\|\mbox{\boldmath{$c$}}\|_{\ell_{2}}},\mbox{ where }\mbox{\boldmath{$\xi$}}\mbox{ is the minimal }\ell_{1}-\mbox{norm}\,\,\mbox{solution of }{\mathcal{A}}\,\mbox{\boldmath{$\xi$}}=\mbox{\boldmath{$c$}}.$$

$$\min\|\mbox{\boldmath{$\rho$}}_{\delta}\|_{\ell_{1}},\hbox{ subject to }{\cal A}\mbox{\boldmath{$\rho$}}_{\delta}=\mbox{\boldmath{$b$}}_{\delta},$$

$${\bf Co}(\mbox{\boldmath{$\rho$}},\mbox{\boldmath{$\rho$}}_{\delta})\leqslant 3\gamma\delta,$$

$${\bf In}(\mbox{\boldmath{$\rho$}},\mbox{\boldmath{$\rho$}}_{\delta})\leqslant 5\gamma\delta.$$

$$\mbox{\boldmath{$\rho$}}_{\ell_{2}}={\cal A}^{*}\mbox{\boldmath{$\mbox{\boldmath$b$}$}},$$

$$\|\mbox{\boldmath$\rho$}\|_{\ell_{1}}\leq\|\mbox{\boldmath{$y$}}\|_{\ell_{1}}<\sqrt{N}\|\mbox{\boldmath{$y$}}\|_{\ell_{2}}<\sqrt{N}\left\|({\mathcal{A}}^{*}_{y}{\mathcal{A}}_{y})^{-1}{\mathcal{A}}^{*}_{y}\right\|_{\ell_{2}}\|\mbox{\boldmath$b$}\|_{\ell_{2}}.$$

$$|\langle\mbox{\boldmath{$a$}}_{i},\mbox{\boldmath{$c$}}_{j}\rangle|<\frac{1}{3M}\,\,\,\,\forall i\hbox{ and }j\,,$$

$$|\langle\mbox{\boldmath{$c$}}_{i},\mbox{\boldmath{$c$}}_{j}\rangle|<\frac{1}{3M}\,\,\,\,\forall i\neq j,$$

$$\gamma ⩽18M^2,$$

$$\forall\,\mbox{\boldmath{$b$}},\,\,\exists\,\,\mbox{\boldmath{$\rho$}}\hbox{ such that }{\cal D}\mbox{\boldmath{$\rho$}}=\mbox{\boldmath{$b$}}\hbox{ and }\|\mbox{\boldmath{$\rho$}}\|_{\ell_{1}}\leqslant\gamma\|\mbox{\boldmath{$b$}}\|_{\ell_{2}}.$$

$$|\langle\mbox{\boldmath{$d$}}_{k},\mbox{\boldmath{$d$}}_{K+s}\rangle|\leqslant\frac{1}{3M},\,\forall k<s+K.$$

$$|\langle\mbox{\boldmath{$d$}}_{k},\mbox{\boldmath{$b$}}\rangle|>\frac{1}{3M}.$$

$$r>\frac{1}{2}\alpha ,where\alpha =\sqrt{1-\frac{1}{9M^2}}.$$

$$\mbox{\boldmath{$b$}}=\sum_{i=1}^{\infty}c_{i}\mbox{\boldmath{$d$}}_{n_{i}}$$

Full text

Keywords: array imaging, $\ell_{1}$-minimization, noise

1 Introduction

In this paper, we are interested in imaging problems formulated as

[TABLE]

so the data vector \mbox{\boldmathb}\in\mathbb{C}^{N} is a linear transformation of the unknown vector \mbox{\boldmath\rho}\in\mathbb{C}^{K} that represents the image. The model matrix ${\cal A}\in\mathbb{C}^{N\times K}$, which is given to us, depends on the geometry of the imaging system and on the sought resolution. Typically, the linear system (1) is underdetermined because only a few linear measurements are gathered, so $N\ll K$. Hence, there exist infinitely many solutions to (1) and, thus, it is a priori not possible to find the correct one without some additional information.

We are interested, however, in imaging problems with sparse scenes. We seek to locate the positions and amplitudes of a small number $M$ of point sources that illuminates a linear array of detectors. This means that the unknown vector $\rho$ is M-sparse, with only a few $M\ll K$ non-zero entries. Under this assumption, (1) falls under the compressive sensing framework [21, 16, 22, 8]. It follows from [16] that the unique M-sparse solution of (1) can be obtained with $\ell_{1}$-norm minimization when the model matrix ${\cal A}$ is incoherent, i.e., when its mutual coherence111The mutual coherence of ${\cal A}$ is defined as \max_{i\neq j}|\langle\mbox{\boldmath{a}}_{i},\mbox{\boldmath{a}}_{j}\rangle|, where the column vectors \mbox{\boldmath{a}}_{i}\in\mathbb{C}^{N} of ${\cal A}$ are normalized to one, so that \|\mbox{\boldmath{a}}_{i}\|_{\ell_{2}}=1 $\forall\,i=1,\ldots,K$. is smaller than $1/(2M)$. The same result can be obtained assuming ${\cal A}$ obeys the M-restricted isometry property [8], which basically states that all sets of M-columns of ${\cal A}$ behave approximately as an orthonormal system.

In our imaging problems these incoherence conditions can be satisfied only for coarse image discretizations that imply poor resolution. To retain resolution and recover the position of the sources with higher precision we propose to extend the theory so as to allow for some coherence in ${\cal A}$. To this end, we show that uniqueness for the minimal $\ell_{1}$-norm solution of (1) can be obtained under less restrictive conditions on the model matrix ${\cal A}$. More specifically, given the columns of ${\cal A}$ that correspond to the support of $\rho$, we define their vicinities as the sets of columns that are almost parallel 222The vicinity of a column \mbox{\boldmath{a}}_{i} is defined as the set of all columns \mbox{\boldmath{a}}_{j} such that |\langle\mbox{\boldmath{a}}_{i},\mbox{\boldmath{a}}_{j}\rangle|\geq 1/(3M). to them. With this definition, our first result set out in Proposition 1 states that if the sources are located far enough from each other, so that their vicinities do not overlap, we can recover their positions exactly with noise free data. Furthermore, in the presence of small noise, their position is still approximately recoverable, in the sense that most of the solution vector is supported in the vicinities while some small noise (grass) is present away from them.

This result finds interesting applications in imaging. As we explain in Section 2, in array imaging we seek to find the position of point sources that are represented as the non-zero entries of $\rho$. Our result states under what conditions the location of these objects can be determined with high precision. It can be also used to explain super-resolution, i.e., the significantly superior resolution that $\ell_{1}$-norm minimization provides compared to the conventional resolution of the imaging system, i.e., the Rayleigh resolution. For instance, super-resolution have been studied using sparsity promotion for sparse spike trains recovery from band-limited measurements. Donoho [15] showed that spike locations and their weights can be exactly recovered for a cutoff frequency $f_{c}$ if the minimum spacing $\Delta$ between spikes is large enough, so $\Delta>1/f_{c}$. Candès and Fernandez-Granda [11] showed that $\ell_{1}$-norm minimization guarantees the exact recovery if $\Delta>1/2f_{c}$. Super-resolution has also been studied for highly coherent model matrices ${\cal A}$ that arise in imaging under the assumption of well separated objects when the resolution is below the Rayleigh threshold [19, 3, 4]. These works include results regarding the robustness of super-resolution in the presence of noise.

Our theory also addresses the robustness to noise of the minimal $\ell_{1}$-norm solution. Specifically, we show that for noisy data the solution $\rho$ can be separated into two parts: (1) the coherent part which is supported inside the vicinities, and (2) the incoherent part, usually referred to as grass, that is small and it is present everywhere. A key observation of our work is that the $\ell_{1}$-images get worse as $\sqrt{N}$ when there is noise in the data and, thus, $\ell_{1}$-norm minimization fails when the number of measurements $N$ is large. This basically follows from (13) in Proposition 1 which relates the $\ell_{1}$ norm of the solution to the $\ell_{2}$ norm of the data, so

[TABLE]

The key quantity here is the constant $\gamma$, which for usual imaging matrices ${\mathcal{A}}$ is proportional to $\sqrt{N}$.

To overcome this problem we introduce in Proposition 2 the noise collector matrix $\mathcal{C}\in\mathbb{C}^{N\times\Sigma}$ and propose to solve instead the augmented linear system [{\cal A}\,|\,{\cal C}]\mbox{\boldmath\rho}=\mbox{\boldmathb}. The dimension of the unknown vector $\rho$ is, thus, augmented by $\Sigma$ components which do not have any physical meaning. They correspond to a fictitious sources that allows us to better approximate the noisy data. The natural question is how to build the noise collector matrix. Theoretically, the answer is given in the proof of Proposition 2 in Section 3, which is constructive. The key is that the column vectors of $[{\cal A}\,|\,{\cal C}]$ form now a frame in which the noisy vector $b$ can be well approximated. As a consequence, we obtain a bound on the constant $\gamma$ ($\gamma<18M^{2}$) which is now independent of $N$. The drawback of this construction is that we need exponentially many vectors, that is $\Sigma\lesssim e^{N}$. This would suggest that the noise collector may not be practical. However, the numerical experiments show that with a large enough number of columns in $\mathcal{C}$ selected at random (as i.i.d. Gaussian random variables with mean zero and variance $1/N$) the $\ell_{1}$-norm minimization problem is regularized and the minimal $\ell_{1}$-norm solution is found.

The paper is organized as follows. In Section 2, we formulate the array imaging problem. In Section 3, we present in a abstract linear algebra framework the conditions under which $\ell_{1}$-minimization provides the exact solution to problem (1) with and without noise. This section contains our main results. In Section 4, we illustrate with numerical simulations how our abstract theoretical results are relevant in imaging sparse sources with noisy data. Section 5 contains our conclusions.

2 Passive array imaging

We consider point sources located inside a region of interest called the image window IW. The goal of array imaging is to determine their positions and amplitudes using measurements obtained on an array of receivers. The array of size $a$ has $N$ receivers separated by a distance $h$ located at positions \vec{\mbox{\boldmath{x}}}_{r}, $r=1,\ldots,N$ (see Figure 1). They can measure single or multifrequency signals with frequencies $\omega_{l}$, $l=1,\dots,S$. The $M$ point sources, whose positions {\vec{\mbox{\boldmath{z}}}}_{j} and complex-valued amplitudes $\alpha_{j}\in\mathbb{C}$, $j=1,\dots,M$, we seek to determine, are at a distance $L$ from the array. The ambient medium between the array and the sources can be homogeneous or inhomogeneous.

In order to form the images we discretize the IW using a uniform grid of points \vec{\mbox{\boldmath{y}}}_{k}, $k=1,\ldots,K$, and we introduce the true source vector

[TABLE]

such that

[TABLE]

We will not assume that the sources lie on the grid, i.e., typically {\vec{\mbox{\boldmath{z}}}}_{j}\neq\vec{\mbox{\boldmath{y}}}_{k} for all $j$ and $k$. To write the data received on the array in a compact form, we define the Green’s function vector

[TABLE]

at location \vec{\mbox{\boldmath{y}}} in the IW, where G(\vec{\mbox{\boldmath{x}}},\vec{\mbox{\boldmath{y}}};\omega) denotes the free-space Green’s function of the homogeneous medium. This function characterizes the propagation of a signal of angular frequency $\omega$ from point \vec{\mbox{\boldmath{y}}} to point \vec{\mbox{\boldmath{x}}}, so (2) represents the signal received at the array due to a point source of amplitude one, phase zero, and frequency $\omega$ at \vec{\mbox{\boldmath{y}}}. If the medium is homogeneous

[TABLE]

The signal received at \vec{\mbox{\boldmath{x}}}_{r} at frequency $\omega_{l}$ is given by

[TABLE]

If we normalize the columns of ${\cal A}$ to one and stack the data in a column vector

[TABLE]

then the source vector $\rho$ solves the system {\cal A}\,\mbox{\boldmath\rho}=\mbox{\boldmathb}, with the $(N\cdot S)\times K$ matrix

[TABLE]

The system {\cal A}\,\mbox{\boldmath\rho}=\mbox{\boldmathb} relates the unknown vector \mbox{\boldmath\rho}\in\mathbb{C}^{K} to the data vector \mbox{\boldmathb}\in\mathbb{C}^{(N\cdot S)}. This system of linear equations can be solved by appropriate $\ell_{2}$ and $\ell_{1}$ methods.

Remark 1**.**

For simplicity of the presentation, we restricted ourselves to the passive array imaging problem where we seek to determine a distribution of sources. The active array imaging problem can be cast under the same linear algebra framework assuming the linearized Born approximation for scattering [12]. In that case, we still obtain a system of the form {\cal A}_{s}\,\mbox{\boldmath\rho}=\mbox{\boldmathb}, where $\rho$ is the reflectivity of the scatterers, $b$ is the data, and ${\cal A}_{s}$ is a model matrix for the scattering problem defined in a similar manner to (6). Even more, when multiple scattering is not negligible the problem can also be cast as in (1); see [13] for details. Therefore, the theory presented in the next sections can be applied to the scattering problems provided that the matrix ${\cal A}_{s}$ satisfies the assumptions of Propositions 1 and 2.

3 $\ell_{1}$ minimization-based methods

In the imaging problems considered here we assume that the sources occupy only a small fraction of the image window IW. This means that the true source vector $\rho$ is sparse, so the number of its entries that are different than zero, denoted by $M$, is much smaller than its length $K$. Thus, we assume M=|\mathrm{supp}(\mbox{\boldmath\rho})|\ll K. This prior knowledge changes the imaging problem substantially because we can exploit the sparsity of $\rho$ by formulating it as an optimization problem which seeks the sparsest vector in $\mathbb{C}^{K}$ that equates model and data. Thus, for a single measurement vector $b$ we solve

[TABLE]

Above, and in the sequel, we denote by $\|\cdot\|_{\ell_{2}}$, $\|\cdot\|_{\ell_{1}}$ and $\|\cdot\|_{\ell_{\infty}}$ the $\ell_{2}$, $\ell_{1}$ and $\ell_{\infty}$ norms of a vector, respectively.

The $\ell_{1}$ minimization problem (7) can be solved efficiently in practice. Several methods have been proposed in the literature to solve it. Here are some of them: orthogonal matching pursuit [7], homotopy [29, 28, 18], interior-point methods [1, 33], gradient projection [20], sub-gradient descent methods in primal and dual spaces [25, 6], and proximal gradient in combination with iterative shrinkage-thresholding [26, 27, 2]. In this work we use GeLMA [24], a semi-implicit version of the primal-dual method [14] that converges to the solution of the following problem independently of the regularization parameter $\tau$: define the function

[TABLE]

for \mbox{\boldmath{\rho}}\in\mathcal{C}^{K} and \mbox{\boldmath{z}}\in\mathcal{C}^{N}, and determine the solution as

[TABLE]

In the literature of compressive sensing we find the following theoretical justification of the $\ell_{1}$-norm minimization approach. If we assume decoherence of the columns of ${\cal A}$, so

[TABLE]

then the $M$-sparse solution of {\cal A}\mbox{\boldmath{\rho}}=\mbox{\boldmath{b}} is unique, and it can be found as the solution of (7) [21, 16, 22]. Numerically, the $\ell_{1}$-norm minimization approach works under less restrictive conditions than the decoherence condition (9) suggests. In fact, our imaging matrices almost never satisfy (9).

Consider a typical imaging regime with central wavelength $\lambda_{0}$. Assume we use $S=36$ equally spaced frequencies covering a bandwidth that is $10\%$ of the central frequency. The size of the array is $a$ and the distance between the array and the IW is $L=a$. An IW of size is $30\lambda_{0}\times 30\lambda_{0}$ is discretized using a uniform grid with mesh size $\lambda_{0}/2\times\lambda_{0}/2$. For such parameters, every column vector \mbox{\boldmath{a}}_{i} has at least sixty two other column vectors \mbox{\boldmath{a}}_{j} so that |\langle\mbox{\boldmath{a}}_{i},\mbox{\boldmath{a}}_{j}\rangle|\geqslant 1/16. Thus, our matrices are fairly far from satisfying the decoherence condition (9) if we want to recover, say, 8 sources. Numerically, however, the $\ell_{1}$ minimization works flawlessly.

Physically, a pair of columns \mbox{\boldmath{a}}_{i} and \mbox{\boldmath{a}}_{j} are coherent, so |\langle\mbox{\boldmath{a}}_{i},\mbox{\boldmath{a}}_{j}\rangle|\approx 1, if the corresponding grid-points in the image are close to each other. In other words, when \mbox{\boldmath{a}}_{i} lies in a vicinity of \mbox{\boldmath{a}}_{j} (and vice versa). We assume, though, that the sources are far apart and, thus, the the set of columns indexed by the support of the true source vector $\rho$ does satisfy the the decoherence condition (9). The above observation motivates the following natural conjecture. Perhaps, the $\ell_{1}$ minimization works well because it suffices to satisfy (9) only on the support of $\rho$. Our main result supports this conjecture.

3.1 Main results

When data is perturbed by small noise, the following qualitative description of the image could be observed. Firstly, some pixels close to the points where the sources are located become visible. Secondly, a few pixels away from the sources are also visible. The latter is usually referred as grass. In order to quantify the observed results we need to modify the the decoherence condition (9) and introduce the vicinities.

Definition 1**.**

Let \mbox{\boldmath{\rho}}\in\mathbb{C}^{K} be an $M$-sparse solution of {\cal A}\mbox{\boldmath{\rho}}=\mbox{\boldmath{b}}, with support $T=\{i:\rho_{i}\neq 0\}$ 333Below and in the rest of the paper the notation $\rho_{i}$ means the $i$th entry of the vector $\rho$. In contrast, we use the notation \mbox{\boldmath{\rho}}_{i} to represent the $i$th vector of a set of vectors.. For any $j\in T$ define the corresponding vicinity of \mbox{\boldmath{a}}_{j} as

[TABLE]

For any vector \mbox{\boldmath{\eta}}\in\mathbb{C}^{K} its coherent misfit to $\rho$ is

[TABLE]

whereas its incoherent remainder with respect to $\rho$ is

[TABLE]

Proposition 1**.**

Let $\rho$ be an $M$-sparse solution of {\cal A}\mbox{\boldmath{\rho}}=\mbox{\boldmath{b}}, and let $T$ be its support. Suppose the vicinities $S_{j}$ from Definition 1 do not overlap, and let $\gamma>0$ be defined as

[TABLE]

Let \mbox{\boldmath{\rho}}_{\delta} be the minimal $\ell_{1}$-norm solution of the noisy problem

[TABLE]

with \|\mbox{\boldmath{b}}-\mbox{\boldmath{b}}_{\delta}\|_{\ell_{2}}\leqslant\delta. Then,

[TABLE]

and

[TABLE]

If $\delta=0$, and $\Upsilon$ does not contain collinear vectors, we have exact recovery: \mbox{\boldmath{\rho}}_{\delta}=\mbox{\boldmath{\rho}}.

Proposition 1 is proved in A. As it follows from this proof, our pessimistic bound $1/(3M)$ could be sharpened to the usual bound (9) found in the literature. We did not strive to obtain sharper results because it will make the proofs more technical and, more importantly, because the concept of vicinities describes well the observed phenomena in imaging with this bound.

When there is no noise so $\delta=0$, Proposition 1 tells us that the M-sparse solution of {\cal A}\mbox{\boldmath{\rho}}=\mbox{\boldmath{b}} can be recovered exactly by solving the $\ell_{1}$ minimization problem under a less stringent condition than (9). Note that we allow for the columns of ${\cal A}$ to be close to collinear. When there is noise so $\delta\neq 0$, this Proposition shows that if the data $b$ is not exact but it is known up to some bounded vector, the solution \mbox{\boldmath{\rho}}_{\delta} of the minimization problem (14) is close to the solution of the original (noiseless) problem in the following sense. The solution \mbox{\boldmath{\rho}}_{\delta} can be separated into two parts: the coherent part supported in the vicinities $S_{j}$ of the true solution, $j\in T$, and the incoherent part, which is small for low noise, and that is supported away from these vicinities. Other stability results can be found in [8, 9, 17, 30, 19, 3].

Let us now make some comments regarding the relevance of this result in imaging. Vicinities, as defined in (10), are related to the classical $\ell_{2}$-norm resolution theory. Indeed, recall Kirchhoff migration imaging given by the $\ell_{2}$-norm solution

[TABLE]

where ${\cal A}^{*}$ is the conjugate transpose of ${\cal A}$. Since \mbox{\boldmathb}={\cal A}\mbox{\boldmath\rho} the resolution analysis of KM relies on studying the behaviour of the inner products |\langle\mbox{\boldmath{a}}_{i},\mbox{\boldmath{a}}_{k}\rangle|. We know from classical resolution analysis [5] that the inner products |\langle\mbox{\boldmath{a}}_{i},\mbox{\boldmath{a}}_{k}\rangle| are large for points \mbox{\boldmath{y}}_{k} that fall inside the support of the KM point spread function, whose size is $\lambda L/a$ in cross-range (parallel to the array) and $c/B$ in range (perpendicular to the array). Hence, we expect the size of the vicinities to be proportional to these classical resolution limits, with an appropriate scaling factor that is inversely proportional to the sparsity $M$. This is illustrated with numerical simulations in Section 4 (right column of Figure 2).

Under this perspective, one could argue that Proposition 1 tells us the well known result that a good reconstruction can be obtained for well separated sources. Proposition 1, however, gives us more information, it provides an $\ell_{1}$-norm resolution theory for imaging: when vicinities do not overlap, there is a single non-zero element of the source associated within each vicinity. Permitting the columns of ${\cal A}$ to be almost collinear inside the vicinities allows for a fine discretization inside the vicinities and therefore the source can be recovered with very high precision. Furthermore, recovery is exact for noiseless data.

The assumptions in Proposition 1 are sufficient conditions but not necessary. Our numerical simulations illustrate exact recovery in more challenging situations, where the vicinities are not well separated (central images in the second row of Figure 2).

For noisy data, Proposition 1 says that it is the concept of vicinities that provide the adequate framework to look at the error between the true solution and the one provided by the $\ell_{1}$-norm minimization approach. Specifically, the error is controlled by the coherent misfit (11) and the incoherent remainder (12), which are shown to be small when the noise is small in $\ell_{2}$. This means that the reconstructed source is supported mainly in the vicinities $S_{j}$ of the true solution, $j\in T$, and the grass in the image is low, i.e., the part of the solution supported away from the vicinities $S_{j}$ is small.

Proposition 1 implies that a key to control the noise is the constant $\gamma$ defined in (13). In general, we have $\gamma=O(\sqrt{N})$. Indeed, let $y$ be the minimum $\ell_{2}$-norm solution of the problem {\cal A}\mbox{\boldmath{\rho}}=\mbox{\boldmath{b}} such that its support has at most size $N$. Let ${\mathcal{A}}_{y}$ be the submatrix of ${\mathcal{A}}$ that contains the columns that correspond to the non-zero entries of $y$. Then, the minimum $\ell_{1}$ solution $\rho$ satisfies (by Cauchy-Schwartz \|\mbox{\boldmath{x}}\|_{\ell_{1}}\leqslant\sqrt{N}\|\mbox{\boldmath{x}}\|_{\ell_{2}},\forall\mbox{\boldmath{x}}\in\mathbb{C}^{N})

[TABLE]

This means that the quality of the image deteriorates as the number of measurements $N\to\infty$. The remedy that we propose to this is to augment the imaging matrix ${\cal A}$ with a “noise collector” ${\cal C}$ as described in the following Proposition.

Proposition 2**.**

There exists a $N\times\Sigma$ noise collector matrix ${\cal C}$, with $\Sigma\lesssim e^{N}$, such that the columns of the augmented matrix ${\cal D}=[{\cal A}\,|\,{\cal C}]$ satisfy \|\mbox{\boldmath{d}}_{j}\|=1,

[TABLE]

[TABLE]

*and there is a positive constant *

[TABLE]

such that

[TABLE]

Proof.

Let \mbox{\boldmath{d}}_{i}=\mbox{\boldmath{a}}_{i}, for $i=1,\dots,K$. We will construct iteratively a sequence of vectors \mbox{\boldmath{d}}_{K+1}=\mbox{\boldmath{c}}_{1}, \mbox{\boldmath{d}}_{K+2}=\mbox{\boldmath{c}}_{2}, $\dots$, \mbox{\boldmath{d}}_{K+\Sigma}=\mbox{\boldmath{c}}_{\Sigma} such that for each $s=1\dots\Sigma$

[TABLE]

The iteration will terminate at a finite step, say, $\Sigma$. At the termination step we will have that for any $b$, \|\mbox{\boldmath{b}}\|_{\ell_{2}}=1 there exists $k\leqslant\Sigma+K$ such that

[TABLE]

The finite time termination is a consequence of a volume growth estimate. Namely, if (19) holds for all $i\neq j\leqslant\Sigma$, then the points \mbox{\boldmath{c}}_{i}, $i=1,2,\dots\Sigma$ are centers of non-overlapping balls of radius $r$. The radius is bounded below:

[TABLE]

Thus the iteration will terminate at a finite step. Furthermore, if $r<\sqrt{2}$ then the number $\Sigma\lesssim e^{N\log{\frac{\sqrt{2}}{r}}}$ as the dimension $N\to\infty$, because $(r/\sqrt{2})^{N}\Sigma\sim 1$.

Let us finally estimate $\gamma$ in (20). Without loss of generality, we may assume \|\mbox{\boldmath{b}}\|_{\ell_{2}}=1. By our construction, there exists $k\leqslant\Sigma+K$ such that (22) holds. Thus we can choose \mbox{\boldmath{d}}_{n_{1}} and $c_{1}$ so that $|c_{1}|\leqslant 1$ and \mbox{\boldmath{b}}_{1}=\mbox{\boldmath{b}}-c_{1}\mbox{\boldmath{d}}_{n_{1}} satisfies \|\mbox{\boldmath{b}}_{1}\|_{\ell_{2}}\leqslant\alpha. Using (22) inductively we can find a sequence \{\mbox{\boldmath{d}}_{n_{i}}\}_{i=1}^{\infty}, and a sequence $\{c_{i}\}_{i=1}^{\infty}$, so that $|c_{i}|\leqslant\alpha^{i-1}$ and the vectors \mbox{\boldmath{b}}_{n}=\mbox{\boldmath{b}}-\sum_{i=1}^{n}c_{i}\mbox{\boldmath{d}}_{n_{i}} satisfy \|\mbox{\boldmath{b}}_{n}\|_{\ell_{2}}\leqslant\alpha^{n}. Therefore,

[TABLE]

and

[TABLE]

by the triangle inequality. ∎

Proposition 2 is an important result as it shows that the constant $\gamma$ in (13) can be made independent of $N$ by augmenting the columns of the linear system with columns of a noise collector matrix ${\cal C}$. The columns of ${\cal C}$ are required to be decoherent to the columns of ${\mathcal{A}}$ (see (18)), and decoherent between them (see (19)). Recalling that the columns of ${\mathcal{A}}$ for the imaging problem are Green’s vectors corresponding to points in the imaging window, we stress that the columns of $\mathcal{C}$ do not admit a physical interpretation. They do not correspond to any points in the imaging window or elsewhere. Similarly, the $\Sigma$ last components of the augmented unknown vector $\rho$ in (21) do not have a physical meaning. They correspond to fictitious auxiliary unknowns that are introduced to regularize the $\ell_{1}$-norm minimization problem.

The drawback in this theory is that the size of the noise collector is exponential $\Sigma\lesssim e^{N}$. This makes it impractical. Our numerical experiments, however, indicate great improvement in the performance of $\ell_{1}$-norm minization with $\Sigma\lesssim 10K$ when the columns of $\mathcal{C}$ are selected at random (its entries are i.i.d. Gaussian random variables with mean zero and variance $1/N$). This works well for additive mean zero uncorrelated noise. For other types of noise, the idea is to construct a library that represents the values that the noise vector \delta\mbox{\boldmathb} takes. It is the elements of this library that should be used as columns of the noise collector matrix $\mathcal{C}$. A different approach can be followed when the noise \delta\mbox{\boldmath{b}} is sparse so its $\ell_{1}$-norm is small. Then, ${\cal C}$ could be simply taken as the $N\times N$ identity matrix $I$. This approach has been proposed and analyzed in [23] and provides exact recovery for sparse noise vectors \delta\mbox{\boldmathb}.

In the next section we present numerical results to illustrate the relevance of our theory in imaging sparse sources. We focus our attention in the case of additive mean zero uncorrelated noise which is not sparse. The results show a dramatic improvement using the noise collector.

4 Imaging results in the framework of propositions 1 and 2

We illustrate here the relevance of Propositions 1 and 2 in imaging. We compare the solution \mbox{\boldmath{\rho}}_{\ell_{1}} obtained with the $\ell_{1}$-norm minimization algorithm GeLMA [24], and the $\ell_{2}$-norm Kirchhoff migration solution (17). Our results illustrate:

The well-known super-resolution for $\ell_{1}$, meaning that \mbox{\boldmath\rho}_{\ell_{1}} determines the support of the unknown $\rho$ with higher accuracy than the conventional resolution limits, provided the assumptions of Proposition 1 are satisfied. 2. 2.

The equally well known sensitivity of $\ell_{1}$ to additive noise. This is made more precise in the imaging context where the constant $\gamma$ in (13) grows with the number of measurements as $\sqrt{NS}$, where $NS$ is the total number of measurements acquired by $N$ receivers at $S$ frequencies. We observe that, for a given level of noise, the $\ell_{1}$-norm reconstruction deteriorates as the number of measurements increases. 3. 3.

The noise collector matrix $\mathcal{C}$ stabilizes $\ell_{1}$-norm minimization in the presence of noise.

We also show how the bandwidth, the array size, and the number of sources affect the vicinities defined in (10). The numerical results are not specialized to a particular physical regime. They illustrate only the role of the Propositions 1 and 2 in solving the associated linear systems.

Imaging setup

The images are obtained in a homogeneous medium with an active array of $N=25$ transducers. We collect measurements corresponding to $S=25$ frequencies equispaced in the bandwidth. Thus, the length of the data vector $b$ is $NS=625$. The ratios between the array size $a$ and the distance $L$ to IW, and between the bandwidth $2B$ and the central frequency $\omega_{0}$ vary in the numerical experiments, so the classical Rayleigh resolution limits change. The size of the IW is fixed. It is discretized using a uniform grid of $K=3721$ points of size $\lambda_{0}/2$ in range and cross-range directions.

The images have been formed by solving the $\ell_{1}$-norm minimization problem (7) using the algorithm GeLMA in [24]. GeLMA is an iterative shrinkage-thresholding algorithm that provides the exact solution, for noiseless data, independently of the value of the parameter $\tau$ (see (8)) used to promote the sparsity of the images.

Results for noiseless data. Super-resolution and $\ell_{1}$-reconstructions

Figure 2 shows the results obtained for a relatively large array and a relatively large bandwidth corresponding to ratios $a/L=1/2$ and $(2B)/\omega_{0}=1/2$ when the data is noiseless. From left to right we show the \mbox{\boldmath{\rho}}_{\ell_{2}} solution (17), the \mbox{\boldmath{\rho}}_{\ell_{1}} solution obtained from (7), the comparison between \mbox{\boldmath{\rho}}_{\ell_{1}} (red stars) and the true solution $\rho$ (green circles), and the vicinities $S_{j}$ defined in (10) plotted with different colors. The top and bottom rows show images with $M=4$ and $M=8$ sources, respectively. The exact locations of the sources are indicated with white crosses in the two leftmost columns. The $M=4$ sources in the top row are very far apart: their vicinities do not overlap as it can be seen in the top right image. In this case, all the conditions of Proposition 1 are satisfied and we find the exact source distribution by $\ell_{1}$-norm minimization. The $M=8$ sources in the bottom row are closer, and their vicinities are larger; according to (10) the size of the vicinities increases with $M$. In fact, their vicinities overlap as it can be seen in the bottom right image. Still, the $\ell_{1}$-norm minimization algorithm finds the exact solution.

The classical resolution limits for this setup are $c_{0}/(2B)=2\lambda_{0}$ in range and $\lambda_{0}L/a=2\lambda_{0}$ in cross-range. This means that the resolution of the $\ell_{2}$-norm solutions is of the order $O(2\lambda_{0})$; see the left column of Figure 2. Recall that our discretization is $\lambda_{0}/2$, that is four times finer than the classical resolution limit. Thus, each source roughly corresponds to a four by four pixel square, which is what the \mbox{\boldmath{\rho}}_{\ell_{2}} solutions show. Note that for $M=8$, because two sources are quite close, the \mbox{\boldmath{\rho}}_{\ell_{2}} solution only displays $7$ sources. The ability of $\ell_{1}$-norm minimization to determine the location of the sources with a better accuracy than the classical resolution limits is referred to as super-resolution.

We stress that if the IW is discretized using a very fine grid, with grid size smaller than the classical resolution limit, then the columns of the matrix ${\cal A}$ are almost parallel and the decoherence condition (9) is violated. The columns that are almost parallel to those indexed by the support of the true solution are contained in the vicinities (10). The number of columns that belong to the vicinities depend on the imaging system. To illustrate the effect of the array and bandwidth sizes on the size of the vicinities we plot in Figure 3 the vicinity of one source for $M=4$. From left to right we plot the vicinities for $[a/L,2B/\omega_{0}]=[1/2,1/2]$, $[a/L,2B/\omega_{0}]=[1/2,1/4]$, $[a/L,2B/\omega_{0}]=[1/4,1/2]$, and $[a/L,2B/\omega_{0}]=[1/4,1/4]$. As expected, the size of the vicinity is proportional to the resolution estimates $\lambda_{0}L/a$ and $c_{0}/(2B)$ in cross-range and range, respectively.

Results for noisy data. Stabilization of $\ell_{1}$-norm minimization using the noise collector matrix ${\cal C}$

We add now mean zero uncorrelated noise to the data. We examine the results for different values of the signal to noise ratio (SNR). We consider first the same imaging configuration as in Figure 2 with $M=4$ sources. The number of data is $NS=625$ and the number of unknowns is $K=3721$. In the top row of Figure 4 we plot the minimal $\ell_{1}$-norm image obtained by solving problem (7) when the SNR is $4$dB. The true solution is shown with white crosses. It is apparent that, even for this moderate level of noise, $\ell_{1}$-norm minimization fails to give a good image.

The problem can be alleviated using the noise collector matrix ${\cal C}$, as it can be seen in the results shown in the bottom row of Figure 4. To construct the noise collector matrix ${\cal C}$ that verifies the assumptions of Proposition 2, we take its columns \mbox{\boldmath{c}}_{j} to be random vectors in $\mathbb{C}^{NS}$ with mean zero and variance $1/(NS)$. Their $\ell_{2}$-norm tends to one as $NS\rightarrow\infty$, and we check that conditions (18) and (19) are satisfied. In theory, the number of columns $\Sigma$ should be very large, of the order of $e^{NS}$, but in practice, we obtain stable results with $\Sigma$ of the order of $10^{4}$, which is roughly $3K$.

The solution \mbox{\boldmath{\mbox{\boldmath$\rho$}}}_{\ell_{1}}\in\mathbb{C}^{K+\Sigma} obtained with the noise collector can be decomposed into two vectors; the vector \mbox{\boldmath{\mbox{\boldmath$\rho$}}}_{iw}\in\mathbb{C}^{K} corresponding to the sought solution in the $IW$, and the vector \mbox{\boldmath{\mbox{\boldmath$\rho$}}}_{noise}\in\mathbb{C}^{\Sigma} that absorbs the noise. We display these two vectors in the bottom right plot of Figure 4. The first $K$ components correspond to \mbox{\boldmath{\mbox{\boldmath$\rho$}}}_{iw} and the remaining $\Sigma$ components to \mbox{\boldmath{\mbox{\boldmath$\rho$}}}_{noise}. It is remarkable that the vector \mbox{\boldmath{\mbox{\boldmath$\rho$}}}_{iw} is very close to the true solution and that it contains only some small grass. This means that both the coherent misfit (15) and the incoherent remainder (16) are now small. This is in accordance with the theoretical error estimates (15) and (16), where $\gamma$ is now independent of the dimension of the data vector $NS$; see (20).

In the next figures, we consider an imaging setup with a large aperture $a/L=1$ and a large bandwidth $(2B)/\omega_{0}=1$. Moreover, we increase the pixel size to $\lambda_{0}$ in both range and cross-range directions, so the Rayleigh resolution is of the order of a pixel. With this imaging configuration, the columns of the model matrix ${\cal A}$ are less coherent than in the previous numerical experiments. We plot in Figure 5 the $\ell_{1}$-norm image for a SNR $=4$dB. With a less coherent matrix ${\cal A}$ the results are very good. This highlights the inherent difficulty in imaging when high resolution is required as in Figure 4 .

For the particular low imaging resolution configuration considered in Figure 5 we obtain good results for a large noise level corresponding to SNR $=0$dB; see the top row of Figure 6 where $N\,S=625$ as before. However, when we increase the number of measurements to $NS=1369$, the image obtained with $\ell_{1}$-norm minimization turns out to be useless; see the bottom row of Figure 6 . This illustrates the counter-intuitive fact that $\ell_{1}$-norm minimization does not always benefit from more data, at less if the data is highly contaminated with noise. This is so because the constant $\gamma$ in (13) depends on the length of the data vector $b$ as $\sqrt{NS}$.

As before, this problem can be fixed with the noise collector as it can be seen in Figure 7. Again, the noise is effectively absorbed for both $NS=625$ (top row) and $NS=1369$ (bottom row) measurements using a matrix collector with a relatively small number of columns, many less than $e^{NS}$ as Proposition 2 suggests.

We finish with one last example that shows that the use of the noise collector makes $\ell_{1}$-norm minimization competitive for imaging sparse scenes because it provides stable results with super-resolution even for highly corrupted data. We consider the example with $M=8$ sources and SNR $=0$dB. The array and the bandwidth are relatively large ($a/L=1/2$, $(2B)/\omega_{0}=1/2$), so the classical $\ell_{2}$-norm resolution is of the order $O(2\lambda_{0})$, as in Figure (4). In Figure 8 we show, from left to right, (i) the minimal $\ell_{1}$-norm solution without noise collector, which fails to give a good image, (ii) the $\ell_{2}$-norm solution (17), which is stable to additive noise but does not resolve nearby sources, and (iii) the minimal $\ell_{1}$-norm solution with the noise collector, which provides a very precise and stable image.

5 Discussion

In this paper, we consider imaging problems that can be formulated as underdetermined linear systems of the form {\cal A}\,\rho_{\delta}=\mbox{\boldmathb}_{\delta}, where ${\cal A}$ is an $N\times K$ model matrix with $N\ll K$, and \mbox{\boldmathb}_{\delta} is the $N$-dimensional data vector contaminated with noise. We assume that the solution is an M-sparse vector in $\mathbb{C}^{K}$, corresponding to the $K$ pixels of the IW. We consider additive noise in the data, so the data vector can be decomposed as \mbox{\boldmathb}_{\delta}=\mbox{\boldmathb}+\delta\mbox{\boldmathb}, where $b$ is the data vector in the absence of noise and \delta\mbox{\boldmathb} is the noise vector. We provide a theoretical framework that allows us to examine under what conditions the $\ell_{1}$-minimization problem admits a solution that is close to the exact one. We also shown that, for our imaging problems, $\ell_{1}$-minimization fails when the noise level is high and the dimension $N$ of the data vector \mbox{\boldmathb}_{\delta} increases. The reason is that the error is proportional to the square root of $N$.

To alleviate this problem and increase the robustness of $\ell_{1}$-minimization, we propose a regularization strategy. In particular, we seek the solution of [{\cal A}\,|\,{\cal C}]\,\mbox{\boldmath\rho}_{\delta}=\mbox{\boldmathb}_{\delta}, where the $N\times\Sigma$ matrix ${\cal C}$ is a noise collector. Thus, the unknown $\rho_{\delta}$ is now a vector in $\mathbb{C}^{K+\Sigma}$. The first $K$ components of the unknown correspond to the distribution of sources in the IW, while the $\Sigma$ next components do not correspond to any physical quantity. They are introduced to provide a fictitious source distribution given by an appropriate linear combination of the columns of ${\cal C}$ that produces a good approximation to \delta\mbox{\boldmathb}. The main idea is to create a library of noises. The columns of the noise collector matrix are elements of this library and they are constructed so as to be incoherent with respect to the columns of ${\cal A}$. Theoretically, the dimension $\Sigma$ of the noise collector increases exponentially with $N$, which suggests that it may not be useful in practice. Our numerical results show, however, robustness for $\ell_{1}$-minimization in the presence of noise when a large enough number of columns $\Sigma\precsim 10K$ is used to build the noise collector matrix.

Our first findings on the noise collector are very encouraging. We have shown that its use improves dramatically the robustness of $\ell_{1}$-norm reconstructions when the data are corrupted with additive uncorrelated noise. Many other questions ought to be addressed. Some directions of our future research concern the following aspects: what happens with other types of noise?, can we design noise collectors adaptively depending on the noise in the data?, what if the noise comes from wave propagation in a random medium?, can we design a noise collector for this case?, how much do we need to know about the noise so as to design a good noise collector?, can we retrieve this information from the data? Some of these questions will be addressed somewhere else.

Acknowledgments

Part of this material is based upon work supported by the National Science Foundation under Grant No. DMS-1439786 while the authors were in residence at the Institute for Computational and Experimental Research in Mathematics (ICERM) in Providence, RI, during the Fall 2017 semester. The work of M. Moscoso was partially supported by Spanish MICINN grant FIS2016-77892-R. The work of A.Novikov was partially supported by NSF grants DMS-1515187, DMS-1813943. The work of C. Tsogka was partially supported by AFOSR FA9550-17-1-0238.