Sparse Recovery Beyond Compressed Sensing: Separable Nonlinear Inverse Problems

TL;DR

This paper develops a theoretical framework for sparse recovery in separable nonlinear inverse problems, extending convex programming methods beyond traditional compressed sensing assumptions, with applications in medical imaging and geophysics.

Contribution

It introduces a new theory for sparse recovery in deterministic measurement settings, proving convex programming can successfully recover parameters under certain correlation conditions.

Findings

Convex programming reliably recovers parameters when they are sufficiently distinct.

The theory applies to deterministic measurement operators violating incoherence conditions.

Numerical experiments demonstrate effectiveness in heat-source localization and brain activity estimation.

Abstract

Extracting information from nonlinear measurements is a fundamental challenge in data analysis. In this work, we consider separable inverse problems, where the data are modeled as a linear combination of functions that depend nonlinearly on certain parameters of interest. These parameters may represent neuronal activity in a human brain, frequencies of electromagnetic waves, fluorescent probes in a cell, or magnetic relaxation times of biological tissues. Separable nonlinear inverse problems can be reformulated as underdetermined sparse-recovery problems, and solved using convex programming. This approach has had empirical success in a variety of domains, from geophysics to medical imaging, but lacks a theoretical justification. In particular, compressed-sensing theory does not apply, because the measurement operators are deterministic and violate incoherence conditions such as the restricted-isometry property. Our main contribution is a theory for sparse recovery adapted to deterministic settings. We show that convex programming succeeds in recovering the parameters of interest, as long as their values are sufficiently distinct with respect to the correlation structure of the measurement operator. The theoretical results are illustrated through numerical experiments for two applications: heat-source localization and estimation of brain activity from electroencephalography data.

Equations

$$f(t):=\sum _{i=1}^k{c}_i{\phi }_t({\theta }_i),$$

$$y$$

$$\underset{{\tilde{\theta }}_1,\dots ,{\tilde{\theta }}_k\in {\mathbb{R}}^p,\tilde{c}\in {\mathbb{R}}^k}{\mathrm{minimize}}{\Vert y-\sum _{i=1}^k{\tilde{c}}_i\vec{\phi }({\tilde{\theta }}_i)\Vert }_2^2.$$

$$\mu :=\sum _{i=1}^k{c}_i{\delta }_{{\theta }_i},$$

$$y={\int }_{{\mathbb{R}}^p}\vec{\phi }(\theta )\mu (\text{d}\theta ).$$

$$\underset{\tilde{\mu }}{\mathrm{minimize}}$$

$${\int }_{{\mathbb{R}}^p}\vec{\phi }(\theta )\tilde{\mu }(\text{d}\theta )=y,$$

$${\Vert \mu \Vert }_{\mathrm{TV}}:=\underset{f\in \mathcal{C},{\Vert f\Vert }_{\infty }\le 1}{\sup }\left|\int f\text{d}\mu \right|,$$

$$\underset{\tilde{\mu }}{\mathrm{minimize}}$$

$${\int }_{{\mathbb{R}}^p}\vec{\phi }(\theta )\tilde{\mu }(\text{d}\theta )=y.$$

$$y$$

$$=\Phi x.$$

$${\rho }_{\theta }(\eta ):=⟨\vec{\phi }(\theta ),\vec{\phi }(\eta )⟩,$$

$${\rho }_{\theta }^{(q,r)}(\eta ):={\vec{\phi }}^{(q)}(\theta {)}^T{\vec{\phi }}^{(r)}(\eta ),$$

$${\rho }_{\theta }^{(0,2)}(\eta )$$

$$|{\rho }_{\theta }^{(1,2)}(\eta )|$$

$$|{\rho }_{\theta }^{(0,0)}(\eta )|$$

$$|{\rho }_{\theta }^{(1,0)}(\eta )|$$

$$|{\rho }_{\theta }^{(0,r)}(\eta )|$$

$$|{\rho }_{\theta }^{(1,r)}(\eta )|$$

$$\underset{i\ne j,{\theta }_i,{\theta }_j\in \Theta }{\min }\left|{\theta }_i-{\theta }_j\right|>2D+\Delta \sigma ,\Delta :=\max \left(\log (1+2(C_{0,0}+C_{1,1})),{\lambda }_1,{\lambda }_2\right),$$

$${\lambda }_1=2\log \left(\frac{2(2(2C_{0,0}+C_{1,1}-C_{1,1}{\gamma }_2+C_{0,1}{\gamma }_3)+1-{\gamma }_2)}{-C_{0,0}+\sqrt{C_{0,0}^2+4(1-{\gamma }_2)(2(2C_{0,0}+C_{1,1}-C_{1,1}{\gamma }_2+C_{0,1}{\gamma }_3)+1-{\gamma }_2)}}\right),$$

$${\lambda }_2=2\log \left(\frac{2(2((2C_{0,0}+C_{1,1}){\gamma }_0+C_{0,2}+C_{0,1}{\gamma }_1)+{\gamma }_0)}{-C_{0,2}+\sqrt{C_{0,2}^2+4{\gamma }_0(2((2C_{0,0}+C_{1,1}){\gamma }_0+C_{0,2}+C_{0,1}{\gamma }_1)+{\gamma }_0)}}\right),$$

$$C_{0,1}{C}_{1,0}$$

$$C_{0,1}{C}_{1,2}$$

$$\gamma ={\left[\begin{array}{cccc}0.185& 0.983& 0.788& 0.868\end{array}\right]}^T,C=\left[\begin{array}{ccc}2.818& 3.348& 4.200\\ 6.786& 8.060& 10.113\end{array}\right],$$

$$d({\theta }_i,{\theta }_j)=\frac{|{\theta }_i-{\theta }_j|-D_i-D_j}{\max ({\sigma }_i,{\sigma }_j)}.$$

$$\Delta >\max \left(\log (1+2(C_{0,0}+C_{1,1})),{\lambda }_1,{\lambda }_2\right),$$

$$y:=\sum _{i=1}^k{c}_i\vec{\phi }({\theta }_i)+z.$$

$$\underset{\tilde{\mu }}{\mathrm{minimize}}$$

$${\Vert {\int }_{{\mathbb{R}}^p}\vec{\phi }(\theta )\tilde{\mu }(\text{d}\theta )-y\Vert }_2\le \xi ,$$

$$\Phi :=\left[\begin{array}{ccc}\vec{\phi }({\eta }_1)& \cdots & \vec{\phi }({\eta }_m)\end{array}\right].$$

$$\underset{\tilde{x}}{\mathrm{minimize}}$$

$$\mathrm{subject}\mathrm{to}$$

$$x_j:=\{\begin{array}{ll}{c}_i& \text{if ηj=θi,}\\ 0& \text{otherwise,}\end{array}$$

Full text

Keywords. Sparse recovery, convex programming, incoherence, correlated measurements, dual certificates, nonlinear inverse problems, source localization.

1 Introduction

1.1 Separable Nonlinear Inverse Problems

The inverse problem of extracting information from nonlinear measurements is a fundamental challenge in many applied domains, including geophysics, microscopy, astronomy, medical imaging, and signal processing. In this work, we focus on separable nonlinear (SNL) problems [44, 45], where the data are modeled as samples from a linear combination of functions that depend nonlinearly on certain quantities of interest. Depending on the application, these quantities may represent neuronal activity in a human brain, oscillation frequencies of electromagnetic waves, locations of fluorescent probes in a cell, magnetic-resonance relaxation times of biological tissues, or positions of celestial objects in the sky.

Mathematically, the goal in an SNL problem is to estimate $k$ parameters $\theta_{1}$, …, $\theta_{k}\in\mathbb{R}^{p}$ from samples of a function

[TABLE]

where $c_{1}$, …, $c_{k}\in\mathbb{R}$ are unknown coefficients. The dependence between each component and the corresponding parameter at a particular value of $t$ is governed by a nonlinear map $\varphi_{t}:\mathbb{R}^{p}\rightarrow\mathbb{R}$. For simplicity of exposition, we assume that $t$ is one-dimensional and $f$ is a real-valued function, but the framework can be directly extended to multidimensional and complex-valued measurements. The data are samples of $f$ at $n$ locations $s_{1}$, …, $s_{n}\in\mathbb{R}$

[TABLE]

where $\vec{\varphi}(\theta_{i})_{j}:=\varphi_{s_{j}}(\theta_{i})$, $1\leq j\leq n$. Each $\vec{\varphi}(\theta_{i})\in\mathbb{R}^{n}$ is a feature vector associated to one of the parameters $\theta_{i}$. Without loss of generality, we assume that the feature vectors are normalized, i.e. $\|\vec{\varphi}(\theta)\|_{2}=1$ for all $\theta\in\mathbb{R}^{p}$. The following examples illustrate the importance of SNL problems in a range of applications.

• •

Deconvolution of point sources: Deconvolution consists of estimating a signal from samples of its convolution with a fixed kernel $K$. When the signal is modeled as a superposition of point sources or spikes, representing fluorescent probes in microscopy [5, 52], celestial bodies in astronomy [75] or interfaces between geological layers in seismography [82], this is an SNL problem where $\theta_{1}$, …, $\theta_{k}$ are the locations of the spikes. In that case, $\varphi_{t}\left(\theta\right)$ is a shifted copy of the convolution kernel $K\left(t-\theta\right)$, as illustrated in the top row of Figure 1.

• •

Spectral super-resolution: Super-resolving the spectrum of a multisinusoidal signal from samples taken over a short period of time is an important problem in communications, radar, and signal processing [85]. This is an SNL problem where $\varphi_{t}\left(\theta\right)$ is a complex exponential $\exp\left(-i2\pi\theta t\right)$ with frequency $\theta$ (see the second row of Figure 1).

• •

Heat-source localization: Finding the position of several heat sources in a material with known conductivity from temperature measurements is an SNL problem where $\varphi_{t}\left(\theta\right)$ is the Green’s function of the heat equation parametrized by the location $\theta$ of a particular heat source [57]. The bottom row of Figure 1 shows an example (see Section 4.1 for more details).

• •

Estimation of neural activity: Electroencephalography measurements of the electric potential field on the surface of the head can be used to detect regions of focalized activity in the brain [68]. The data are well approximated by an SNL model where the parameters are the locations of these regions [64]. The function $\varphi_{t}\left(\theta\right)$ represents the potential at a specific location $t$ on the scalp, which originates from neural activity at position $\theta$ in the brain. This function can be computed by solving the Poisson differential equation taking into account the geometry and electric properties of the head [70]. Figure 2 shows an example. See Section 4.2 for more details.

• •

Quantitative magnetic-resonance imaging: The magnetic-resonance relaxation times $T_{1}$ and $T_{2}$ of biological tissues govern the local fluctuations of the magnetic field measured by MR imaging systems [69]. MR fingerprinting is a technique to estimate these parameters by fitting an SNL model where each component corresponds to a different tissue [60, 63, 90]. In this case, the parameter $\theta\in\mathbb{R}^{2}$ encodes the values of $T_{1}$ and $T_{2}$ and the function $\varphi_{t}\left(\theta\right)$ can be computed by solving the Bloch differential equations [7].

1.2 Reformulation as a Sparse-Recovery Problem

A natural approach to estimate the parameters of an SNL model is to solve the nonlinear least-squares problem,

[TABLE]

Unfortunately, the resulting cost function is typically nonconvex and has local minima, as illustrated by the simple example in Figure 3. Consequently, local-descent methods do not necessarily recover the true parameters, even in the absence of noise, and global optimization becomes intractable unless $k$ is very small.

Alternatively, we can reformulate the SNL problem as a sparse-recovery problem and leverage $\ell_{1}$-norm minimization to solve it. This approach was pioneered in the 1970s by geophysicists working on spike deconvolution in the context of reflection seismology [92, 23, 55, 80, 25]. Since then, it has been applied to many SNL problems such as acoustic sensing [100, 4], radar [74, 91], electroencephalography (EEG) [83, 98], positron emission tomography (PET) [49, 76, 51], direction of arrival [61, 8], quantitative magnetic resonance imaging  [63, 90], and source localization [57, 62, 72]. Our goal is to provide a theory of sparse recovery via convex optimization explaining the empirical success of this approach.

Let us represent the parameters $\theta_{1}$, …, $\theta_{k}\in\mathbb{R}^{p}$ of an SNL model as a superposition of Dirac measures or spikes in $\mathbb{R}^{p}$, interpreted as a $p$-dimensional parameter space, (1.2),

[TABLE]

where $\delta_{\theta_{i}}$ denotes a Dirac measure supported at $\theta_{i}$. Intuitively, the atomic measure $\mu$ is a signal that encodes the parameters of interest and their corresponding coefficients. The data described by Equation 1.2 can now be expressed as

[TABLE]

The SNL problem is equivalent to recovering $\mu$ from these linear measurements. The price we pay for linearizing is that the linear inverse problem is extremely underdetermined: $y$ has dimension $n$, but $\mu$ lives in a continuous space of infinite dimensionality! To solve the problem, we need to exploit the assumption that the data only depends on a small number of parameters or, equivalently, that $\mu$ is sparse.

For an SNL problem to be well posed, $\theta_{1}$, …, $\theta_{k}$ should be the only set of $k$ or less parameters such that Equation 1.2 holds. In that case, $\mu$ is the solution to the sparse-recovery problem

[TABLE]

where minimization occurs over the set of measures in $\mathbb{R}^{p}$. The cardinality of the support of an atomic measure is a nonconvex function, which is notoriously challenging to minimize. A fundamental insight underlying many approaches to sparse estimation in high-dimensional statistics and signal processing is that one can bypass this difficulty by replacing the nonconvex function with a convex counterpart. In particular, minimizing the $\ell_{1}$ norm instead of the cardinality function has proven to be very effective in many applications.

In order to analyze the application of $\ell_{1}$-norm minimization to SNL problems, we consider a continuous setting, where the optimization variable is a measure supported on a continuous domain. The goal is to obtain an analysis that is valid for arbitrarily fine discretizations of the domain. This is important because, as we will see below, a fine discretization results in a highly-correlated linear operator, which violates the usual assumptions made in the literature on sparse recovery.

In the case of measures supported on a continuous domain, the continuous counterpart of the $\ell_{1}$ norm is the total-variation (TV) norm [79, 40]111Not to be confused with the total variation of a piecewise-constant function used in image processing.. Indeed, the TV norm of the atomic measure $\mu$ in Equation 1.4 equals the $\ell_{1}$ norm of its coefficients $\left\|c\right\|_{1}$. Just as the $\ell_{1}$ norm is the dual norm of the $\ell_{\infty}$ norm, the TV norm is defined by

[TABLE]

where the supremum is taken over all continuous functions in the unit $\mathcal{L}_{\infty}$-norm ball. Replacing the cardinality function by this sparsity-promoting norm yields the following convex program

[TABLE]

The goal of this paper is to understand when the solution to Problem (1.8) exactly recovers the parameters of an SNL model.

1.3 Compressed Sensing

Section 1.2 shows that solving an SNL problem is equivalent to recovering a sparse signal from linear underdetermined measurements. This is reminiscent of compressed sensing [27, 20, 42]. In its most basic formulation, the goal of compressed sensing is to estimate a signal $x\in\mathbb{R}^{m}$ with $k$ nonzeros from linear measurements $y\in\mathbb{R}^{n}$ given by $y:=Ax$, where $A\in\mathbb{R}^{n\times m}$ and $m>n$. Remarkably, exact recovery of $x$ is still possible under certain conditions on the matrix $A$, even though the linear system is underdetermined.

Overdetermined linear inverse problems are said to be ill posed when the measurement matrix is ill conditioned. This occurs when there exist vectors that lie close to the null space of the matrix, or equivalently when a subset of its columns is highly correlated. Analogously, the compressed-sensing problem is ill posed if any sparse subset of columns is highly correlated, because this implies that sparse vectors lie close to the null space. Early works on compressed sensing derive recovery guarantees assuming a bound on the maximum correlation between the columns of the measurement matrix $A$ (sometimes called incoherence). They prove that tractable algorithms such as $\ell_{1}$-norm minimization and greedy techniques achieve exact recovery as long as the maximum correlation is of order $n^{-1/2}$ for sparsity levels $k$ of up to order $\sqrt{n}$ [28, 31, 48, 93], even if the data are corrupted by additive noise [30, 94]. These results were subsequently strengthened to sparsity levels of order $n$ (up to logarithmic factors) [27, 67, 18, 6, 12] under stricter assumptions on the conditioning of sparse subsets of columns in the measurement matrix, such as the restricted-isometry property [18] or the restricted-eigenvalue condition [6].

The question is whether compressed-sensing theory applies to SNL problems. Let us consider an SNL problem where the parameter space is discretized to yield a finite-dimensional version of the sparse-recovery problem described in Section 1.2. The measurement model in Equation 1.5 can then be expressed as

[TABLE]

$\Phi$ is a measurement matrix whose columns correspond to the feature vectors of $\eta_{1}$, …, $\eta_{m}\in\mathbb{R}^{p}$, which denote the $m$ points of the discretized parameter space. The signal $x\in\mathbb{R}^{m}$ is the discrete version of $\mu$ in Equation 1.4: a sparse vector, such that $x_{i}=c_{j}$ when $\eta_{j}=\theta_{i}$ for some $i\in\left\{1,\ldots,k\right\}$ and $x_{j}=0$ otherwise. For compressed-sensing theory to apply here, the intercolumn correlations of $\Phi$ should be very low. Figure 4 shows the intercolumn correlations of a typical compressed-sensing matrix: each column is almost uncorrelated with every other column. Figures 5 and 6 show the correlation function $\rho_{\theta}:\mathbb{R}^{p}\rightarrow\mathbb{R}$

[TABLE]

for the different SNL problems discussed in Section 1.1. The intercolumn correlations of the measurement matrix $\Phi$ in the corresponding discretized SNL problems are given by samples of $\rho_{\eta_{i}}$ at the locations of the remaining grid points $\eta_{j}$, $i\neq j$. The contrast with the intercolumn structure of the compressed-sensing matrix is striking: nearby columns in all SNL measurement matrices are very highly correlated. This occurs under any reasonable discretization of the parameter space; discretizing very coarsely would result in inaccurate parameter estimates, defeating the whole point of solving the SNL problem.

The difference between the correlation structure of the measurement matrix in compressed-sensing and SNL problems is not surprising. The entries of compressed-sensing matrices are random. As a result, small subsets of columns are almost uncorrelated with high probability. In contrast, matrices in discretized SNL problems arise from a deterministic model tied to an underlying continuous parameter space and to a function $\varphi_{t}$ that is typically smooth. Since $\varphi_{t}\left(\theta\right)\approx\varphi_{t}\left(\theta^{\prime}\right)$ when $\theta\approx\theta^{\prime}$, nearby columns are highly correlated. These matrices do not satisfy any of the properties of the conditioning of sparse submatrices commonly assumed in compressed sensing. In conclusion, the answer to our previous question is a resounding no: compressed-sensing theory does not apply to SNL problems.

1.4 Beyond Sparsity and Randomness: Separation and Correlation Decay

The fact that compressed-sensing theory does not apply to SNL problems involving deterministic measurements is not a theoretical artifact. Sparsity is not a strong enough condition to ensure that such SNL problems are well posed. If $\varphi_{t}$ is smooth, which is usually the case in applications, the features $\vec{\varphi}\left(\theta\right)$ corresponding to parameters that are clustered in the parameter space are highly correlated. This can be seen in the correlation plots of Figures 5 and 6. As a result, different sparse combinations of features may yield essentially the same data. For a detailed analysis of this issue in the context of super-resolution and deconvolution of point sources we refer the reader to Section 3.2 in [14] (see also [65] and [84]) and Section 2.1 in [3], respectively.

Additional assumptions beyond sparsity are necessary to establish recovery guarantees for SNL problems. At the very least, the features $\vec{\varphi}\left(\theta_{1}\right)$, …, $\vec{\varphi}\left(\theta_{k}\right)$ in the data cannot be too correlated. For arbitrary SNL problems it is challenging to define simple conditions to preclude this from happening. However, in most practical situations, SNL problems exhibit correlation decay, meaning that the correlation function $\rho_{\theta}$ defined in Equation 1.11 is bounded by a decaying function away from $\theta$. This is a natural property: the more separated two parameters $\theta$ and $\theta^{\prime}$ are in the parameter space, the less correlated we expect their features $\vec{\varphi}\left(\theta\right)$ and $\vec{\varphi}\left(\theta^{\prime}\right)$ to be. All the examples in Section 1.1 have correlation decay (see Figures 5 and 6).

For SNL problems with correlation decay there is a simple way of ensuring that the features corresponding to the true parameters $\theta_{1}$, …, $\theta_{k}$ are not highly correlated: imposing a minimum separation between them in the parameter space. The main contribution of this paper is showing that this is in fact sufficient to guarantee that TV-norm minimization achieves exact recovery, under some additional conditions on the derivatives of the correlation function.

1.5 Organization

In Section 2 we propose a theoretical framework for the analysis of sparse estimation in the context of SNL inverse problems. We focus on the case $p=1$ for simplicity, but our results can be extended to higher dimensions, as described in Section 3.5. Our main results are Theorems 2.4 and 2.6, which establish exact-recovery results for SNL problems with correlation decay under a minimum separation on the true parameters. Section 3 contains the proof of these results, which are based on a novel dual-certificate construction. Section 4 illustrates the theoretical results through numerical experiments for two applications: heat-source localization and estimation of brain activity from electroencephalography data.

2 Main Results

2.1 Correlation decay

In this section we formalize the notion of correlation decay by defining several conditions on the correlation function $\rho_{\theta}$ and on its derivatives. Throughout we assume that the problem is one dimensional ($p:=1$).

To alleviate notation we define

[TABLE]

for $q=0,1$ and $r=0,1,2$, where $\vec{\varphi}^{\,(q)}$ is the $q$th derivative of $\vec{\varphi}$. Recall that we assume $\|\vec{\varphi}(\theta)\|_{2}=1$ for all $\theta\in\mathbb{R}$. This implies $\rho_{\theta}(\theta)=1$ and $\rho^{(1,0)}_{\theta}(\theta)=\rho^{(0,1)}_{\theta}(\theta)=0$ for all $\theta\in\mathbb{R}$. Plots of these derivatives are shown in Figure 8 for the deconvolution, super-resolution, and heat-source localization problems. Our conditions take the form of bounds in different regions of the parameter space: a near region, an intermediate region, and a far region, as depicted in Figure 7.

In the near region the correlation can be arbitrarily close to one, but is locally bounded by a quadratic function.

Condition 2.1** (Near Condition).**

The correlation function $\rho_{\theta}$ satisfies the near condition if

[TABLE]

hold for all $\eta$ in $[N^{-},N^{+}]$, where $N^{\pm}:=\theta\pm N$, and $\gamma_{0},\gamma_{1}$ are positive constants.

Equation (2.2) requires correlations to be concave locally, which is natural since the maximum of $\rho_{\theta}$ is attained at $\theta$. Equation (2.3) is a regularity condition that requires $\rho_{\theta}^{(0,2)}(\eta)$ to vary smoothly as we change the center $\theta$. The normalization quantity $\left\|\vec{\varphi}^{\,(1)}(\theta)\right\|_{2}^{2}$ captures how sensitive the features $\vec{\varphi}$ are to perturbations. If this quantity is small for some $i$ then we require more regularity from $\rho_{\theta}$ because $\theta$ is harder to distinguish from nearby points using the measurements $\vec{\varphi}$.

In the intermediate region the correlation function $\rho_{\theta}$ is bounded but can otherwise fluctuate arbitrarily. In addition, we require a similar constraint on its derivative with respect to the position of the center $\theta$.

Condition 2.2** (Intermediate Condition).**

The correlation function $\rho_{\theta}$ satisfies the intermediate condition if

[TABLE]

hold for $\eta<N^{-}$ and $\eta>N^{+}$, where $N^{\pm}$ are defined as in the near condition, and $\gamma_{2},\gamma_{3}$ are positive constants.

In the decay region the correlation and its derivatives are bounded by a decaying function.

Condition 2.3** (Decay Condition).**

The correlation function $\rho_{\theta}$ satisfies the decay condition with decay constant $\sigma>0$ if

[TABLE]

hold for $\eta<D^{-}$ and $\eta>D^{+}$, where $r=0,1,2$, $D^{\pm}:=\theta\pm D$, and $C_{q,r}$ are positive constants.

The choice of exponential decay is for concreteness, and can be replaced by a different summable decay bound222In the case of super-resolution, the decay is not summable, but can be made summable by applying a window to the data, which is standard practice in spectral super-resolution [50].. Figure 8 shows the derivatives of the correlation functions for several SNL problems.

The sample locations $s_{1},\ldots,s_{n}$ have a direct effect on the correlation function, and hence on whether an SNL problem satisfies Conditions 2.1, 2.2 and 2.3. In Figure 9 we depict the correlation function $\rho_{\theta}$ for different sampling patterns in a Gaussian deconvolution problem with $\vec{\varphi}(\theta)_{i}:=\exp(-(\theta-s_{i})^{2}/2)$. For both regular and irregular sampling patterns, the correlation functions are well behaved as long as there are samples close to the true parameter. In contrast, when there are few samples, or the samples are distant from the true parameter, the conditions in Conditions 2.1, 2.2 and 2.3 may be violated. We refer the interested reader to [3] for a more precise analysis of sampling patterns for deconvolution of point sources.

2.2 Exact Recovery for SNL Problems with Uniform Correlation Decay

In this section we focus on SNL problems where the correlation function $\rho_{\theta}$ of the measurement operator is approximately translation invariant, meaning that $\rho_{\theta}$ has similar properties for any value of $\theta$. Examples of such SNL problems include super-resolution, deconvolution, and heat-source localization if the conductivity is approximately uniform. We prove that TV-norm minimization recovers a superposition of Dirac measures exactly as long as the support satisfies a separation condition related to the decay properties of the correlation function and its derivatives.

Theorem 2.4**.**

Let $\Theta:=\{\theta_{1},\ldots,\theta_{k}\}$ be the support of the measure $\mu$ defined in Equation 1.4. Assume that the correlation functions $\rho_{\theta_{i}}$, $\theta_{i}\in\Theta$, satisfy Conditions 2.1, 2.2 and 2.3 for fixed constants $N$, $D$, $\sigma$, $\gamma\in\mathbb{R}^{4}$, and $C\in\mathbb{R}^{2\times 3}$. Then $\mu$ is the unique solution to Problem (1.8) as long as the minimum separation of $\Theta$ satisfies,

[TABLE]

where

[TABLE]

[TABLE]

and the constants $C_{q,r}$ are chosen so that

[TABLE]

Note that the condition in Equation 2.11 is only needed to simplify the statement and proof of our results.

Theorem 2.4 establishes that TV minimization recovers the true parameters of an SNL problem when the support separation is larger than a constant that is proportional to the rate of decay of the correlation function and its derivatives. This separation is measured from the edges of the intermediate regions, as depicted in Figure 10. In stark contrast to compressed-sensing theory, the result holds for correlation functions that are arbitrarily close to one in the near regions, and may have arbitrary bounded fluctuations in the intermediate regions.

Our statement and proof of Theorem 2.4 focuses on clarity over sharpness. For example, in the Gaussian deconvolution problem depicted in Figure 9 with dense uniform samples, the theorem requires a minimum separation of approximately $6.6$ to guarantee exact parameter recovery. To obtain this bound, we let $N=1$, $D=0$, $\sigma=1$,

[TABLE]

and verify the conditions numerically. This is in contrast with the minimum separation of approximately $3$ proven in [3]. Part of this discrepancy comes from applying exponential decay bounds to a Gaussian-shaped correlation function.

Our result requires conditions on the correlation functions centered at the true support $\Theta$, and also on their derivatives. The decay conditions on the derivatives constrain the correlation structure of the measurement operator when we perturb the position of the true parameters. For example, they implicitly bound pairwise correlations centered in a small neighborhood of the support. Exploring to what extent these conditions are necessary is an interesting question for future research.

2.3 Exact Recovery for SNL Problems with Nonuniform Correlation Decay

The measurement operators associated to many SNL problems of practical interest have nonuniform correlations. Figures 5 and 6 show that this is the case for heat-source localization with spatially-varying conductivity, and for estimation of brain-activity from EEG data. Our goal in this section is to establish recovery guarantees for such problems.

The conditions on the correlation structure of the measurement operator required by Theorem 2.4 only pertain to the correlation functions centered at the true parameters. In order to generalize the result, we allow the correlation function centered at each parameter to satisfy the conditions in Section 2.1 with different constants. This makes it possible for the correlation to have near and intermediate regions of varying widths around each element of the support, as well as different decay constants in the far region. Our main result is that TV-norm minimization achieves exact recovery for SNL problems with nonuniform correlation structure, as long as the support satisfies a minimum-separation condition dependent on the corresponding support-centered correlation functions.

Let $\Theta$ be the support of our signal of interest, and assume $\rho_{\theta_{i}}$ satisfies the decay conditions in Section 2.1 with parameters $\sigma_{i}$, $D_{i}$, and $N_{i}$, which are different for all $\theta_{i}\in\Theta$. Extending the notation, we let $N_{i}^{\pm}:=\theta_{i}\pm N_{i}$ and $D_{i}^{\pm}:=\theta_{i}\pm D_{i}$ denote the endpoints of the near and decay regions, respectively. Intuitively, when $\sigma_{i}$ and $D_{i}$ are small, the corresponding correlation function $\rho_{\theta_{i}}$ is “narrower” and should require less separation than “wider” correlation functions with large values of $\sigma_{i}$ and $D_{i}$. This is illustrated in Figure 11, where we depict $\rho_{\theta_{i}}$ for the heat-source localization problem at three different values of $i$. The decay becomes more pronounced towards the right due to the changing thermal conductivity of the underlying medium. For the problem to be well posed, one would expect $\theta_{1}$ to require more separation from other active sources than $\theta_{2}$, which in turn should require more than $\theta_{3}$. We confirm this intuition through numerical experiments in Section 4.1. To make it mathematically precise, we define the following generalized notion of support separation.

Definition 2.5** (Generalized support Separation).**

Suppose for all $\theta_{i}\in\Theta$ that $\rho_{\theta_{i}}$ satisfies Condition 2.3 with parameters $D_{i}$ and $\sigma_{i}$. Define the normalized distance $d(\theta_{i},\theta_{j})$ for $i\neq j$ by

[TABLE]

Assume that $\Theta$ is ordered so that $\theta_{1}<\theta_{2}<\cdots<\theta_{{k}}$. $\Theta$ has separation $\Delta>0$ if $d(\theta_{i},\theta_{j})>|i-j|\Delta$ for all $\theta_{i},\theta_{j}\in\Theta$ with $i\neq j$.

The normalized distance $d(\theta_{i},\theta_{j})$ is measured between the edges of the decay regions of $\theta_{i}$ and $\theta_{j}$, and normalized by the level of decay. This allows sharply decaying correlation functions to be in close proximity with one another. We require $d(\theta_{i},\theta_{j})>|i-j|\Delta$ to prevent the parameters from becoming too clustered. If we only require the weaker condition $d(\theta_{i},\theta_{j})>\Delta$, and if $\sigma_{i}$ grows very quickly with $i$, then we could have $d(\theta_{1},\theta_{j})\approx\Delta$ for all $j>i$. This causes too much overlap between the correlation functions.

Figure 13 gives an example of parameters and correlation functions that satisfy the conditions of Definition 2.5. The following theorem establishes exact-recovery guarantees under a condition on the generalized support separation.

Theorem 2.6**.**

Suppose that for all $\theta_{i}\in\Theta$ that $\rho_{\theta_{i}}$ satisfies Conditions 2.1, 2.2 and 2.3 and (2.11) with constants $N:=N_{i}$, $D:=D_{i}$, $\sigma:=\sigma_{i}$, $C$, and $\gamma$. Note that $C$ and $\gamma$ are the same for each $\theta_{i}$. Then the true measure $\mu$ defined in (1.4) is the unique solution to Problem (1.8) when $\Theta$ has separation $\Delta$ (as determined by Definition 2.5) satisfying

[TABLE]

(2.9), (2.10), and (2.11).

The proof of Theorem 2.6, which implies Theorem 2.4, is given in Section 3. The theorem establishes that TV minimization recovers the true parameters of SNL problems with nonuniform correlation decays when the generalized support separation is larger than a constant. Equivalently, exact recovery is achieved as long as each true parameter $\theta_{i}$ is separated from the rest by a distance that is proportional to the rate of decay of the correlation function centered at $\theta_{i}$. The separation is measured from the edges of the intermediate regions, which can also vary in width as depicted in Figure 10. The result matches our intuition about SNL problems: the parameters can be recovered as long as they yield measurements that are not highly correlated. As mentioned previously, the theorem requires decay conditions on the derivatives of the correlation function, which constrain the correlation structure of the measurement operator.

2.4 Robustness to Noise

In practice, measurements are always corrupted by noisy perturbations. Noise can be taken into account in our measurement model (1.2) by incorporating an additive noise vector $z\in\mathbb{R}^{n}$:

[TABLE]

To adapt the TV-norm minimization problem (1.8) to such measurements, we relax the data consistency constraint from an equality to an inequality:

[TABLE]

where $\xi>0$ is a parameter that must be tuned according to the noise level. Previous works have established robustness guarantees for TV-norm minimization applied to specific SNL problems such as super-resolution [37, 15] and deconvolution [3] at small noise levels. These proofs are based on dual certificates. Combining the arguments in [37, 3] with our dual-certificate construction in Section 3 yields robustness guarantees for general SNL problems in terms of support recovery at high signal-to-noise regimes. We omit the details, since the proof would essentially mimic the ones in [37, 3]. In Figure 14 we show an application of Equation 2.16 to a noisy deconvolution problem taken from [3].

2.5 Discretization

The continuous optimization problem (2.16) can be solved by applying $\ell_{1}$-norm minimization after discretizing the parameter space. This is a very popular approach in practice for a variety of SNL problems [61, 92, 100, 83, 90]. If the true parameters lie on the discretization grid, then our exact-recovery results translate immediately. The following corollary is a discrete version of Theorem 2.6.

Corollary 2.7**.**

Assume $\Theta$ lies on a known discretized grid $G:=\{\eta_{1},\ldots,\eta_{m}\}$ so that $\Theta\subset G$. Furthermore, suppose the conditions of Theorem 2.6 hold so that $\mu$ as defined in (1.4) is the unique solution to Problem (1.8). Define the dictionary $\Phi\in\mathbb{R}^{n\times m}$ by

[TABLE]

Then the $\ell_{1}$-norm minimization problem

[TABLE]

has a unique solution $x\in\mathbb{R}^{m}$ satisfying

[TABLE]

for $j=1,\ldots,m$.

Proof.

If the support of $\tilde{\mu}$ in Problem (1.8) is restricted to lie on $G$, then Problems (1.8) and (2.18) are the same. Thus $\left\|\mu\right\|_{\operatorname{TV}}$ must be smaller than $\left\|\tilde{x}\right\|_{1}$ for any $\tilde{x}$ such that $\Phi\tilde{x}=y$. By assumption $\mu$ is supported on $G$, so the result follows. ∎

Of course, the true parameters may not lie on the grid used to solve the $\ell_{1}$-norm minimization problem. The proof techniques used to derive robustness guarantees for super-resolution and deconvolution in [37, 3] can be leveraged to provide some control over the discretization error. Performing a more accurate analysis of discretization error for SNL problems is an interesting direction for future research.

2.6 Related Work

2.6.1 Sparse Recovery via Convex Programming from Deterministic Measurements

In [29], Donoho and Elad develop a theory of sparse recovery from generic measurements based on the spark and coherence of the measurement matrix $\Phi$. The spark is defined to be the smallest positive value $s$ such that $\Phi$ has $s$ linearly dependent columns. The coherence, which we denote by $M(\Phi)$, is the maximum absolute correlation between any two columns of $\Phi$. The authors show that exact recovery is achieved by $\ell_{1}$-minimization when the number of true parameters is less than $(1+1/M(\Phi))/2$. As discussed in Section 1.3, these arguments are inapplicable to the finely-discretized parameter spaces occurring in SNL problems since neighboring columns of $\Phi$ have correlations approaching one. In [32], the authors provide a support-dependent condition for exact recovery of discrete vectors. Using our notation, they require that $\beta(\Theta)/(1-\alpha(\Theta))<1$ where

[TABLE]

Here $\vec{\varphi}(\eta_{j})$ for $\eta_{j}\notin\Theta$ ranges over the columns of $\Phi$ that do not correspond to the true parameters. This condition is also inapplicable for matrices arising in our problems of interest because $\beta(\Theta)$ approaches one (or larger) for finely-discretized parameter spaces. Sharper exact-recovery guarantees in subsequent works [22, 27, 67, 18, 6, 12] require randomized measurements, and therefore do not hold for deterministic SNL problems as explained in Section 1.3.

2.6.2 Convex Programming Applied to Specific SNL Problems

In [14, 38], the authors establish exact recovery guarantees for super-resolution via convex optimization by leveraging parameter separation (see Section 1.4). Subsequent works build upon these results to study the robustness of this methodology to noise [89, 37, 15, 34], missing data [87], and outliers [39]. A similar analysis is carried out in [3] for deconvolution. The authors establish a sampling theorem for Gaussian and Ricker-wavelet convolution kernels, which characterizes what sampling patterns yield exact recovery under a minimum-separation condition. Other works have analyzed the Gaussian deconvolution problem under nonnegativity constraints [81, 35], and also for randomized measurements [73]. All of these works exploit the properties of specific measurement operators. In contrast, the present paper establishes a general theory that only relies on the correlation structure of the measurements. The works that are closer to this spirit are [2, 88], which analyze deconvolution via convex programming for generic convolution kernels. The results in [2] require quadratically decaying bounds on the first three derivatives of the convolution kernel. In [88], the authors prove exact recovery assuming bounds on the first four derivatives of the autocorrelation function of the convolution kernel. In contrast to these works, our results allow for discrete irregular sampling and for measurement operators that are not convolutional, which is necessary to analyze applications such as heat-source localization or estimation of brain activity.

Algorithms for solving the convex programs arising from SNL problems divide into two categories: grid-based and grid-free (or off-the-grid). In grid-based methods the parameter space is discretized, thus yielding a finite-dimensional $\ell_{1}$ minimization problem that can be solved using standard methods (as we have done in Section 4). Grid-free methods attempt to directly solve the infinite-dimensional, but with weaker guarantees (see [53] and [58] for review articles). Recently, several authors [26, 36, 10] have proposed grid-free algorithms based on the conditional gradient method [43].

2.6.3 Other Methodologies

SNL parameter recovery can be formulated as a nonlinear least squares problem [45]. For a fixed value of the parameters $\theta_{1},\ldots,\theta_{k}$, the optimal coefficients $c_{1},\ldots,c_{k}$ in (1.1) have a closed form solution. This makes it possible to minimize the nonlinear cost function with respect to $\theta_{1},\ldots,\theta_{k}$ directly, a technique known as variable projection [44]. As shown in Figure 3, a downside to this approach is that it may converge to suboptimal local minima, even in the absence of noise.

Prony’s method [24, 85] and the finite-rate-of-innovation (FRI) framework [97, 33] can be applied to tackle SNL problems, as long one can recast them as spectral super-resolution problems. This provides a recovery method that avoids discretizing the parameter space. The FRI framework has also been applied to arbitrary non-bandlimited convolution kernels [95] and nonuniform sampling patterns [71], but without exact-recovery guarantees. These techniques have been recently extended by Dragotti and Murray-Bruce [66] to physics-driven SNL problems. By approximating complex exponentials with weighted sums of Green’s functions, they are able to recast parameter recovery as a related spectral super-resolution problem that approximates the true SNL problem.

3 Proof of Exact-Recovery Results

3.1 Dual Certificates

To prove Theorem 2.6 we construct a certificate that guarantees exact recovery.

Proposition 3.1** (Proof in Appendix A).**

Let $\Theta=\{\theta_{1},\ldots,\theta_{k}\}\subset\mathbb{R}$ denote the support of the signed atomic measure $\mu$. Assume that for any sign pattern $\xi\in\{\pm 1\}^{k}$ there is a $\tilde{q}\in\mathbb{R}^{p}$ such that $\widetilde{Q}(\theta):=\tilde{q}^{T}\vec{\varphi}(\theta)$ satisfies

[TABLE]

Then $\mu$ is the unique solution to problem (1.8).

To prove exact recovery of a signal we need to show that it is possible to interpolate the sign pattern of its amplitudes, which we denote by $\xi$, on its support $\Theta$ using an interpolating function $\widetilde{Q}(\theta)$ that is expressible as a linear combination of the coordinates of $\vec{\varphi}(\theta)$. The coefficient vector of this linear combination, denoted by $\tilde{q}$, is known as a dual certificate in the literature because it certifies recovery and is a solution to the Lagrange dual of problem (1.8):

[TABLE]

Dual certificates have been widely used to derive guarantees for inverse problems involving random measurements, including compressed sensing [19, 17], matrix completion [11] and phase retrieval [13]. In such cases, the certificates are usually constructed by leveraging concentration bounds and other tools from probability theory [96]. In contrast, our setting is completely deterministic. More recently, dual certificates have been proposed for specific deterministic problems such as super-resolution [14] and deconvolution [3]. Our goal here is to provide a construction that is valid for generic SNL models with correlation decay.

3.2 Correlation-Based Dual Certificates

Our main technical contribution is a certificate that only depends on the correlation function of the measurement operator. In contrast, certificate constructions in previous works on SNL problems [14, 3, 2] typically rely on problem-specific structure, with the exception of [88] which proposes a certificate for time-invariant problems based on autocorrelation functions.

In our SNL problems of interest the function $\vec{\varphi}(\theta)$ mapping the parameters of interest to the data is assumed to be continuous and smooth. As a result, Equations 3.1 and 3.2 imply that any valid interpolating function $\widetilde{Q}$ reaches a local extremum at each $\theta_{i}\in\Theta$. Equivalently, $\widetilde{Q}$ satisfies the following $2k$ interpolation equations for all $\theta_{i}\in\Theta$:

[TABLE]

Inspired by this observation, we define

[TABLE]

where $\alpha_{i},\beta_{i}\in\mathbb{R}$, $i=1,\ldots,k$, are chosen so that $Q(\theta):=q^{T}\vec{\varphi}(\theta)$ satisfies the interpolation equations. Crucially, this choice of coefficients yields an interpolation function $Q$ that is a linear combination of correlation functions centered at $\Theta$,

[TABLE]

In essence, we interpolate the sign pattern of the signal on its support using support-centered correlations. Since $\rho_{\theta_{i}}(\theta_{i})=1$ and $\rho_{\theta_{i}}^{(1,0)}(\theta_{i})=0$, Equations 3.4 and 3.8 imply $\alpha_{i}\approx\xi_{i}$ when $\rho(\theta,\eta)$ and $\rho^{(1,0)}(\theta,\eta)$ decay as $|\theta-\eta|$ grows large and the support is sufficiently separated. The term that depends on $\beta$ can be interpreted as a correction to the derivatives of $Q$ so that (3.5) is satisfied. The normalizing factor used in (3.8) makes this explicit:

[TABLE]

for $i=1,\ldots,k$. Figure 15 illustrates the construction for the heat-source localization problem. The construction is inspired by previous interpolation-based certificates tailored to super-resolution [14] and deconvolution [3].

In the remainder of this section we show that our proposed construction yields a valid certificate if the conditions of Theorem 2.6 hold. In Section 3.3 we prove Lemma 3.2, which establishes that the interpolation equations have a unique solution and therefore $Q$ satisfies (3.1).

Lemma 3.2**.**

Under the assumptions of Theorem 2.6 there exist $\alpha,\beta\in\mathbb{R}^{k}$ which uniquely solve Equations 3.4 and 3.5.

In Section 3.4 we prove Lemma 3.3, which establishes that $Q$ satisfies (3.2). This completes the proof of Theorem 2.6.

Lemma 3.3**.**

Let $\alpha,\beta\in\mathbb{R}^{k}$ be the coefficients obtained in Lemma 3.2. Under the assumptions of Theorem 2.6 , the interpolating function $Q$ defined in (3.8) satisfies $|Q(\theta)|<1$ for $\theta\notin\Theta$.

3.3 Proof of Lemma 3.2

To simplify notation we define the $i$th normalized correlation and its derivatives by

[TABLE]

for $q=0,1$ and $r=0,1,2$, where $\rho_{\theta_{i}}^{(q,r)}$ is defined in Equation 2.1. Using this notation we have

[TABLE]

To prove Lemma 3.2 we establish the following stronger result, which also gives useful bounds on $\alpha$ and $\beta$. Throughout we assume that $\rho$ satisfies Conditions 2.1, 2.2 and 2.3 and (2.11) with parameters $\gamma$, $C$, $N$, $D$, and $\sigma$, and that $\Theta$ has separation $\Delta$.

Lemma 3.4**.**

Let $s:=2e^{-\Delta}/(1-e^{-\Delta})$. If $\Delta>\log(1+2(C_{0,0}+C_{1,1}))$ then there are $\alpha,\beta\in\mathbb{R}^{k}$ which uniquely solve equations (3.4) and (3.5). Furthermore, we have

[TABLE]

To prove Lemma 3.4 we begin by rewriting Equations (3.4) and (3.5) in block matrix-vector form for the $Q$ function defined in Equation 3.8:

[TABLE]

where $\mathcal{I}\in\mathbb{R}^{{k}\times{k}}$ is the identity matrix, and $P^{(q,r)}\in\mathbb{R}^{{k}\times{k}}$ satisfies

[TABLE]

for $i\neq j$ and $P^{(q,r)}_{ii}=0$. To see why $P_{ii}^{(1,0)}=0$ agrees with equations (3.4) and (3.5), note that

[TABLE]

The atoms are normalized– $\left\|\vec{\varphi}(\theta)\right\|_{2}=1$ for all $\theta$– which implies

[TABLE]

For the same reason it follows that $\hat{\rho}_{\theta_{i}}^{(0,1)}(\theta_{i})=0$.

Our plan is to bound the norm of each $P^{(q,r)}$. If these norms are sufficiently small then the matrix in Equation 3.13 is nearly the identity, and the desired result follows from a linear-algebraic argument. Define $\epsilon_{i}^{(q,r)}(\theta)$ by

[TABLE]

where $q=0,1$ and $r=0,1,2$. Here we think of $\theta$ as a point close to $\theta_{i}$, so $\epsilon_{i}^{(q,r)}(\theta)$ captures the cumulative correlation from the other, more distant elements of $\Theta$. We expect $\epsilon_{i}^{(q,r)}(\theta)$ to be small when $\rho$ has sufficient decay and $\Theta$ is well separated. For a matrix $A$ let $\left\|A\right\|_{\infty}$ denote the infinity-norm defined by

[TABLE]

$\left\|A\right\|_{\infty}$ equals the maximum sum of absolute values in any row of $A$. We have

[TABLE]

The following lemma shows that equation (3.13) is invertible when $\epsilon^{(q,r)}$ is sufficiently bounded.

Lemma 3.5** (Proof in Appendix B).**

Suppose

[TABLE]

Then the matrix in (3.13) is invertible and

[TABLE]

To apply Lemma 3.5 we must first bound $\epsilon^{(q,r)}$ for $q,r\in\{0,1\}$. By the decay and separation conditions (Condition 2.3 and Definition 2.5), $\theta_{j}$ lies in the exponentially decaying tail of $\rho_{\theta_{i}}$, for $i\neq j$. This gives

[TABLE]

where $s:=2e^{-\Delta}/(1-e^{-\Delta})$ and the distance $d$ is defined in Definition 2.5. As this bound is independent of $i$, we have

[TABLE]

for $q,r\in\{0,1\}$. In terms of the conditions of Lemma 3.5, we obtain

[TABLE]

assuming Equation 2.11 ($C_{1,0}C_{0,1}=C_{0,0}C_{1,1}$) holds. Thus we have $c<1$, as required by (3.20), when $s<\frac{1}{C_{0,0}+C_{1,1}}$. Using this, we can find conditions on $\Delta$ so that the hypotheses of Lemma 3.5 hold. If $\Delta>\log(1+\kappa)$ for some $\kappa>0$ then, using that $f(x)=e^{-x}/(1-e^{-x})$ is decreasing for $x>0$, we have

[TABLE]

Letting $\kappa=2C_{1,1}$ this shows that $\Delta>\log(1+2C_{1,1})$ implies

[TABLE]

giving (3.19). For condition (3.20) we let $\kappa=2(C_{0,0}+C_{1,1})$. Then (3.30) shows that $\Delta>\log(1+2(C_{0,0}+C_{1,1}))$ implies

[TABLE]

as required. Thus when $\Delta>\log(1+2(C_{0,0}+C_{1,1}))$ both conditions of Lemma 3.5 hold. By plugging (3.28) and (3.29) into the bounds of Lemma 3.5 we obtain Lemma 3.4. This completes the proof of Lemma 3.2.

3.4 Proof of Lemma 3.3

Lemma 3.2 implies that we can solve (3.4) and (3.5) for $\alpha$ and $\beta$ and then obtain $Q$ via (3.8). To prove Lemma 3.3 we must show that $|Q(\theta)|<1$ for $\theta\in\Theta^{c}$. This is accomplished in two steps. We first show that $|Q(\theta)|<1$ for $\theta$ that are not in the near region of any $\rho_{\theta_{i}}$, $i=1,\ldots,k$. Secondly, we show that for any $\theta$ in the near region of some $\rho_{\theta_{i}}$ we have $Q^{(2)}(\theta)<0$, where we assume $\xi_{i}=1$ without loss of generality. This proves that $Q$ has a local maximum at $\theta_{i}$, and is smaller than one nearby.

To bound $Q$ and $Q^{(2)}$ we apply the triangle inequality to (3.8), obtaining the following lemma.

Lemma 3.6**.**

Fix $\theta_{i}\in\Theta$ and let $Q$ be defined as in (3.8). For all $\theta\in\mathbb{R}$

[TABLE]

If $\hat{\rho}_{\theta_{i}}^{(0,2)}(\theta)\leq 0$ and $\xi\in\{-1,+1\}^{k}$ is the sign pattern interpolated by $Q$ then

[TABLE]

In the next lemma we show that the separation conditions in Definition 2.5 ensure that the support does not cluster together. We assume that $\theta_{1}<\theta_{2}<\cdots<\theta_{|\Theta|}$.

Lemma 3.7** (Proof in Appendix C).**

Fix $\theta_{i}\in\Theta$ and let $\theta>\theta_{i}$. If $i\leq k-1$, assume that

[TABLE]

Then

[TABLE]

where $\Delta$ is defined in Definition 2.5, as long as Conditions 2.1, 2.2 and 2.3 and (2.11) hold.

Inequality (3.35) implies that $\theta_{i}$ is the closest element of $\Theta$ to $\theta$, if we use the generalized distance normalized by $\sigma$.

Bounding $|Q(\theta)|$ Outside the Near Region

Our goal is to bound $|Q(\theta)|$ for $\theta\in\mathbb{R}$ that are not in the near region of any $\rho_{\theta_{i}}$. We can assume, flipping the axis if necessary, that the conditions of Lemma 3.7 hold for the $\theta_{i}$ closest to $\theta$ and that $\theta\geq N_{i}^{+}$ (recall that $N_{i}^{+}=\theta_{i}+N_{i}$ is the boundary of the near region of $\rho_{\theta_{i}}$). Then $\theta$ lies in the exponentially decaying tail of $\rho_{\theta_{j}}$ for $j\neq i$, so that

[TABLE]

where $x=e^{-\Delta/2}$ and $s$ is defined in Lemma 3.4. Plugging this into (3.33) and combining with Lemma 3.4 yields

[TABLE]

Solving for $|Q(\theta)|<1$ we obtain

[TABLE]

Isolating $s$ and $x$,

[TABLE]

where we apply Equation 2.11 ($C_{1,0}C_{0,1}=C_{0,0}C_{1,1}$) to cancel terms. Since $s=\frac{2x^{2}}{1-x^{2}}$ we obtain the inequality

[TABLE]

where $a:=2C_{0,0}+C_{1,1}-C_{1,1}\gamma_{2}+C_{0,1}\gamma_{3}>0$, $b:=C_{0,0}$, and $c:=1-\gamma_{2}>0$. The following lemma shows that this inequality is satisfied by our assumptions, completing this part of the proof.

Lemma 3.8** (Proof in Appendix D).**

Let $a,b,c,\Delta>0$, and $x:=e^{-\Delta/2}$. Then the inequality

[TABLE]

is satisfied when

[TABLE]

Bounding $Q^{(2)}(\theta)$ in the Near Region

For the final part of the proof we must prove $Q^{(2)}(\theta)<0$ for $\theta$ in the near region of some $\theta_{i}$ with $\xi_{i}=1$. Since the interpolation equations (3.4) and (3.5) guarantee that $Q(\theta_{i})=1$, strict concavity implies that $Q$ has a unique maximum on $[N_{i}^{-},N_{i}^{+}]$ thus establishing that $Q(\theta)<1$ for $\theta\in[N_{i}^{-},N_{i}^{+}]\setminus\{\theta_{i}\}$. We cannot have $Q(\theta)\leq-1$ on the near region since we already showed that $|Q(\theta)|<1$ outside the near region. We can assume, without loss of generality, that the conditions of Lemma 3.7 hold for some $i$ and that $\theta\leq N_{i}^{+}$. By the same argument used to obtain Equation 3.42, we have

[TABLE]

Plugging this into (3.33) and applying Lemma 3.4 we obtain

[TABLE]

where $\hat{\rho}_{\theta_{i}}^{(0,2)}(\theta)\leq 0$ by Condition 2.1. Solving for $Q^{(2)}(\theta)<0$ yields

[TABLE]

Grouping into terms involving $s$ and $x$ we obtain

[TABLE]

where we apply Equation 2.11 ($C_{0,1}C_{1,2}=C_{0,2}C_{1,1}$) to cancel terms. Letting $a:=(2C_{0,0}+C_{1,1})\gamma_{0}+C_{0,2}+C_{0,1}\gamma_{1}>0$, $b:=C_{0,2}$, and $c:=\gamma_{0}>0$ we obtain an inequality of the form (3.48). Applying Lemma 3.8 completes the proof.

3.5 Extensions to Higher Dimensions

In this section we briefly describe an extension of our dual certificate construction to settings where the parameter space has dimension $p>1$. We leave a more detailed analysis to future work. In higher dimensions, one can build the interpolating function $Q$ by setting

[TABLE]

where $\alpha\in\mathbb{R}^{k}$, $\beta\in\mathbb{R}^{k\times p}$, $\partial_{l}$ denotes the partial derivative with respect to the $l$th coordinate, and $\partial_{1,l}\rho_{\theta}(\eta)$ denotes the partial derivative of $\rho_{\theta}(\eta)$ with respect to the $l$th coordinate of $\theta$. The coefficients are chosen so that $Q$ satisfies the analog of the interpolation Equations 3.4 and 3.5,

[TABLE]

In Figure 16 we show an example of an interpolating function $Q$ for the electroencephalography brain-activity localization problem. The interpolating function is associated to the signal with three active sources of brain activity from Figure 2. To control $Q$ for problems with correlation decay, one can extend the correlation-decay conditions in Section 2.1 to higher dimensions by requiring analogous bounds on the first and second partial derivatives. For example, the local quadratic bound in Equation 2.2 becomes a bound on the eigenvalues of the Hessian of $\rho_{\theta}$. Similar conditions have been used in [14, 2] to obtain recovery guarantees for the super-resolution and deconvolution problems in two dimensions.

4 Numerical Experiments

4.1 Heat-source localization

Our theoretical results establish that convex programming yields exact recovery in parameter-estimation problems that have correlation decay. In this section we investigate this phenomenon numerically for a heat-source localization problem in one dimension. The heat sources are modeled as a collection of point sources,

[TABLE]

where $\Theta$ is a finite number of points in the unit interval and $c_{i}\in\mathbb{R}$ for $i=1,\ldots,k$.

The heat $u(\theta,t)$ at position $\theta$ and time $t$ is assumed to evolve following the heat equation with Neumann boundary conditions (see e.g. [57]),

[TABLE]

on the unit interval, where $c(\theta)$ represents the conductivity of the medium at $\theta$ (see Figure 17). The data are heat measurements $u(j,T)$, where $j$ is sampled on a regular grid of 100 points at a fixed time $T:=10^{-4}$. Our goal is to recover the initial heat distribution at $t=0$ represented by the heat sources $\mu$. This is an SNL problem where $\vec{\varphi}(\theta)$ corresponds to the measurements caused by a source located at $\theta$.

If the conductivity is constant, then $\vec{\varphi}(\theta)$ has a Gaussian-like shape with fixed width sampled on the measurement grid, and the correlation function $\rho_{\theta}$ is also shaped like a Gaussian with fixed width. In this example, the conductivity varies (see Figure 17), which results in correlations that are still mostly Gaussian-like but have very different widths (see Figure 18). The measurement operator therefore has nonuniform correlation decay properties. According to the theory presented in Section 2.3, we expect convex programming to recover the heat-source location as long as they are separated by a minimum separation that takes into account the structure of the correlation decay. To verify whether this is the case we consider two different separation measures. The first just measures the minimum separation between the sources,

[TABLE]

The second takes into account the correlation function of the SNL problem,

[TABLE]

where $\sigma_{i}$, $1\leq i\leq k$, is the standard deviation of the best Gaussian upper bound on the correlation function centered at $\theta_{i}$,

[TABLE]

The question is which of these separations characterizes the performance of the convex-programming approach more accurately. We consider two different heat-source configurations: one where the sources are uniformly spaced and another where they are clustered, as depicted in Figure 18. In the clustered configuration, groups of spikes are placed in regions of low conductivity where each $\rho_{\theta_{i}}$ exhibits sharp decay. This creates a noticeable discrepancy between the two measures of separation. We expect exact recovery to occur for much smaller values of $\Delta_{\operatorname{sep}}$ in the clustered case than it does in the uniform case since it doesn’t account for the correlation decay. In contrast, we expect $\Delta_{\operatorname{corr}}$ to have similar thresholds for exact recovery in both cases, as it captures the inherent difficulty in the problem. To recover these sources we solve the discretized $\ell_{1}$-norm minimization problem described in Section 2.5 using CVX [47]. The measurement matrix is computed by solving the differential equation on a grid of $10^{3}$ points in the parameter space.

Figure 19 shows the results. For both types of support, we observe that exact recovery occurs as long as the minimum separation $\Delta_{\operatorname{sep}}$ between sources is larger than a certain value. However, this value is very different for the two types of support. It is much smaller for the clustered support. Intuitively, this is due to the fact that the clustered sources are mostly located in the region where the conductivity is lower, and consequently the correlation decay is much faster. Our theory suggests that exact recovery will occur at smaller separations for such configurations, which is what we observe. Quantifying separation using $\Delta_{\operatorname{corr}}$ accounts for the variation in correlation decay, resulting in a similar phase transition for both types of support. This is consistent with the theory for SNL problems with nonuniform correlation decay developed in Section 2.3.

4.2 Estimation of Brain Activity via Electroencephalography

In this section we consider the SNL problem of brain-activity localization from electroencephalography (EEG) data. Areas of focalized brain activity are usually known as sources. In EEG, as well as in magnetoencephalography, sources are usually modeled using electric dipoles, represented mathematically as point sources or Dirac measures at the corresponding locations [1]. EEG measurements are samples of the electric potential on the surface of the head obtained using an array of electrodes. The mapping from the source positions and amplitudes to the EEG data is governed by Poisson’s equation. For a fixed model of the head geometry– obtained for example from 3D images acquired using magnetic-resonance imaging– the mapping can be computed by solving the differential equation numerically. By linearity, this corresponds to an SNL model where $\theta$ represents the position of a dipole and $\varphi_{t}\left(\theta\right)$ is the value of the corresponding electrical potential at a location $t$ of the scalp.

To simulate realistic EEG data with associated ground-truth source positions, we use Brainstorm [86], an open-source software for the analysis of EEG and other electrophysiological recordings. We use template ICBM152, which is a nonlinear average of 152 3D magnetic-resonance image scans [41], to model the geometry of the brain and head. The sources are modeled as electrical dipoles situated on the cortical surface, discretized as a tesselated grid with 15,000 points. The orientation of the dipoles is assumed to be perpendicular to the cortex. The data-acquisition sensor array has 256 channels (HydroCel GSN, EGI). The data corresponding to each point on the grid are simulated numerically using the open-source software OpenMEEG [46]. The computation assumes realistic electrical properties for the brain, skull and scalp. As a result, we obtain a $256\times 15,000$ measurement operator that can be used to generate data corresponding to different combinations of sources. Figure 2 shows three examples.

The correlation structure of the EEG linear measurement operator is depicted in Figure 6. The complex geometric structure of the folds in cortex induce an intricate correlation structure close to the source location. However, the correlation presents a clear correlation decay; points that are far enough from the source have low correlation values. Our theoretical results consequently suggest that $\ell_{1}$-norm minimization should succeed in recovering superpositions of sources that are sufficiently separated. This is supported by the fact that for such superpositions, we can build valid dual certificates as illustrated in Figure 16.

In order to test our hypothesis, we randomly choose superpositions of three sources located at approximately the same distance from each other, for a range of distances. For each superposition we generate 8 different measurement vectors by assigning every possible combination of positive and negative unit-norm amplitudes to each source. We then estimate the source locations by solving an $\ell_{1}$-norm minimization problem with equality constraints using CVX [47]. Figure 20 shows the results. When the separation is sufficiently large, exact recovery indeed occurs for all possible patterns of positive and negative amplitudes, as expected from our theoretical analysis. As the separation decreases, recovery fails for some of the patterns, as shown in the image on the right.

5 Conclusion and Future Work

In this work, we establish that deterministic SNL problems can be solved via a tractable convex optimization program, as long as the parameters have a minimum separation in the parameter space with respect to the correlation structure of the measurement operator. As mentioned in Sections 2.4, our results can be used to establish some robustness guarantees at low-noise levels in terms of support recovery. Deriving more precise stability analysis in the spirit of [56] is an interesting research direction. Another interesting open problem is characterizing the performance of reweighting techniques, which are commonly applied in practice to enhance solutions obtained from $\ell_{1}$-norm minimization [21, 90].

A drawback of the sparse-recovery framework discussed in this paper is computational complexity: for SNL problems in two or three dimensions it is very computationally expensive to solve $\ell_{1}$-norm minimization problems, even if we discretize the domain. In the last few years, algorithms that minimize nonconvex cost functions directly via gradient descent have been shown to provable succeed in solving underdetermined inverse problems involving randomized measurements [16, 59]. As illustrated in Figure 3, nonlinear least-squares cost functions associated to deterministic SNL problems often have non-optimal local minima. An intriguing question is how to design cost functions for deterministic SNL problems that can be tackled directly by descent methods.

Another computationally efficient alternative is to perform recovery using machine learning. Recent works suggest calibrating a feedforward network to output the model parameters using a training set of examples, with applications in point-source deconvolution [9] and super-resolution of line spectra [54]. Understanding under what conditions such techniques can be expected to yield accurate estimates is a challenging question for future research.

Acknowledgements

B.B. is supported by the MacCracken Fellowship, and the Isaac Barkey and Ernesto Yhap Fellowship. C.P. is supported by NIH award NEI-R01-EY025673. C.F. is supported by NSF award DMS-1616340.

Appendix A Proof of Proposition 3.1

The following proof is standard, but given here for completeness.

.

Let $\nu$ be feasible for problem (1.8) and define $h=\nu-\mu$. By taking the Lebesgue decomposition of $h$ with respect to $|\mu|$ we can write

[TABLE]

where $h_{\Theta}$ is absolutely continuous with respect to $|\mu|$, and $h_{\Theta^{c}}$, $|\mu|$ are mutually singular. In other words, the support of $h_{\Theta}$ is contained in $\Theta$, and $h_{\Theta^{c}}(\Theta)=0$. This allows us to write

[TABLE]

for some $b\in\mathbb{R}^{|\Theta|}$. Set $\xi:=\operatorname*{sign}(b)$, where we arbitrarily choose $\xi_{j}=\pm 1$ if $b_{j}=0$. By assumption there exists a corresponding $Q$ interpolating $\xi$ on $\Theta$. Since $\mu$ and $\nu$ are feasible for problem (1.8) we have $\int\varphi(\theta)\,dh(\theta)=0$. This implies

[TABLE]

Applying the triangle inequality, we obtain

[TABLE]

where the last inequality is strict if $\left\|h_{\Theta^{c}}\right\|_{\operatorname{TV}}>0$ since $|Q(\theta)|<1$ for $\theta\in\Theta^{c}$. This establishes that $\mu$ is optimal for problem (1.8), and that any other optimal solution must be supported on $\Theta$. Equation A.3 implies that any feasible solution supported on $\Theta$ must be equal to $\mu$ (since $\left\|h_{\Theta}\right\|_{\operatorname{TV}}=0$), completing the proof of uniqueness. ∎

Appendix B Proof of Lemma 3.5

.

For any matrix $A\in\mathbb{R}^{n\times n}$ such that $\left\|A\right\|_{\infty}<1$ the Neumann series $\sum_{j=0}^{\infty}A^{j}$ converges to $(\mathcal{I}-A)^{-1}$, which implies that $\mathcal{I}-A$ is invertible [77]. By the triangle inequality and the submultiplicativity of the $\infty$-norm, this gives

[TABLE]

Setting $A=-P^{(1,1)}$ and applying $\|P^{(1,1)}\|_{\infty}=\epsilon^{(1,1)}<1$ proves $\mathcal{I}+P^{(1,1)}$ is invertible. Let $\mathcal{C}$ be the Schur complement of $\mathcal{I}+P^{(1,1)}$ in (3.13) so that

[TABLE]

By the triangle inequality and (B.1) applied with $A=-P^{(1,1)}$ we obtain

[TABLE]

proving $\mathcal{C}$ is invertible. As both $\mathcal{I}+P^{(1,1)}$ and its Schur complement $\mathcal{C}$ are invertible, the matrix in (3.13) is also invertible (see e.g. [99]), which establishes the first claim.

By applying blockwise Gaussian elimination we solve (3.13) in terms of $\mathcal{C}$ to obtain

[TABLE]

Applying (B.1) and noting that $\left\|\xi\right\|_{\infty}=1$ we obtain the required bounds on $\left\|\alpha\right\|_{\infty}$ and $\left\|\beta\right\|_{\infty}$. Finally,

[TABLE]

which implies $\left\|\alpha-\xi\right\|_{\infty}\leq c\left\|\alpha\right\|_{\infty}$ and completes the proof. ∎

Appendix C Proof of Lemma 3.7

.

Throughout we use the fact that $\Theta$ having separation $\Delta>0$ (Definition 2.5) implies $D_{i}^{+}\leq D_{j}^{-}$ for $j>i$. When $i\leq k-1$, (3.35) implies that

[TABLE]

where the righthand side is the $\sigma$-weighted average of $D_{i}^{+}$ and $D_{i+1}^{-}$. If $\theta\geq D_{i}^{+}$ this gives

[TABLE]

If $\theta<D_{i}^{+}$ then

[TABLE]

by the separation conditions, proving (3.36). If $i+1<j\leq|\Theta|$ then $d(\theta_{i+1},\theta_{j})\geq\Delta(j-(i+1))$, so that

[TABLE]

by the separation conditions applied to $\theta_{j}$ and $\theta_{i+1}$, which implies (3.37). Finally, if $j<i$ then

[TABLE]

by the separation conditions applied to $\theta_{j}$ and $\theta_{i}$, which establishes (3.38). ∎

Appendix D Proof of Lemma 3.8

.

Multiplying through by $1-x^{2}$ (which is positive by assumption) in (3.48) we obtain

[TABLE]

The above inequality is implied by the simpler quadratic inequality

[TABLE]

where the omitted term $-x^{3}b$ is always negative. Since the inequality is satisfied for $x=0$ we obtain the condition

[TABLE]

Translating this to a statement on $\Delta$ we obtain

[TABLE]

∎