Harnessing Sparsity over the Continuum: Atomic Norm Minimization for Super Resolution

TL;DR

This paper explains how atomic norm minimization, a convex optimization method, can be used for super resolution tasks, providing mathematical foundations, performance guarantees, and applications in microscopy.

Contribution

It offers a comprehensive exposition of atomic norm minimization for super resolution, including theoretical guarantees and practical application insights.

Findings

Atomic norm minimization provides effective super resolution solutions.

The approach has proven performance guarantees.

Application demonstrated in fluorescence microscopy image reconstruction.

Abstract

Convex optimization recently emerges as a compelling framework for performing super resolution, garnering significant attention from multiple communities spanning signal processing, applied mathematics, and optimization. This article offers a friendly exposition to atomic norm minimization as a canonical convex approach to solve super resolution problems. The mathematical foundations and performances guarantees of this approach are presented, and its application in super resolution image reconstruction for single-molecule fluorescence microscopy are highlighted.

Equations

$$\mathbf{x}=\sum _i{c}_i{\mathbf{a}}_i,{\mathbf{a}}_i\in \mathcal{A},$$

$${\Vert \mathbf{x}\Vert }_{\mathcal{A}}\triangleq \inf \left\{t\ge 0:\mathbf{x}\in t\cdot \mathrm{conv}\left(\mathcal{A}\right)\right\},$$

$${\Vert \mathbf{x}\Vert }_{\mathcal{A}}=\inf \left\{\sum _i{c}_i:\mathbf{x}=\sum _i{c}_i{\mathbf{a}}_i,c_i>0,{\mathbf{a}}_i\in \mathcal{A}\right\},$$

$$x(t)=\sum _{k=1}^r{c}_k\delta (t-{\tau }_k).$$

$$G(f)=0\text{mbox}whenever|f|>B/2$$

$$y(t)$$

$$Y_{\ell }=G_{\ell }\cdot X_{\ell }+E_{\ell }=G_{\ell }\cdot \left(\sum _{k=1}^r{c}_k{e}^{-j2\pi \ell {\tau }_k}\right)+E_{\ell },$$

$$\mathbf{y}=\text{mbox}diag(\mathbf{g})\mathbf{x}+\mathbf{ϵ},$$

$$\mathbf{z}=\mathbf{x}+\widetilde{\mathbf{ϵ}},$$

$$\mathbf{x}=\sum _{k=1}^r{c}_k\mathbf{a}({\tau }_k),$$

$$\mathbf{a}(\tau )={\left[1,e^{j2\pi \tau },\dots ,e^{j2\pi (n-1)\tau }\right]}^{\top },\tau \in [0,1).$$

$$\mathbf{x}={\int }_0^1\mathbf{a}(t)\mathrm{d}x(t).$$

$${\mathcal{A}}_0=\left\{\mathbf{a}(\tau ):\tau \in \left[0,1\right)\right\},$$

$${\mathcal{A}}_{1D}=\left\{e^{j\varphi }\mathbf{a}(\tau ):\tau \in \left[0,1\right),\varphi \in \left[0,2\pi \right)\right\}.$$

$$\min \Vert x(t){\Vert }_{\mathsf{TV}}\text{mbox}s.t.\mathbf{x}={\int }_0^1\mathbf{a}(t)\mathrm{d}x(t),$$

$$\displaystyle\left\|{\boldsymbol{x}}\right\|_{{\mathcal{A}}}=\inf_{\begin{subarray}{c}{\boldsymbol{u}}\in\mathbb{C}^{n}\\ t>0\end{subarray}}\Big{\{}\frac{1}{2n}\mathop{\rm Tr}$$

$$\displaystyle\begin{bmatrix}\mathop{\rm toep}\left({\boldsymbol{u}}\right)&{\boldsymbol{x}}\\ {\boldsymbol{x}}^{\mathsf{H}}&t\end{bmatrix}\succeq 0\Big{\}}.$$

$$\mathbf{G}\in \mathcal{G}(R)⟺\mathrm{Tr}(\mathbf{G})=\mathrm{Re}(r_0),\sum _{i=1}^{n-k}{G}_{i,i+k}=\frac{r_k}{2},k=1,\dots ,n-1.$$

$$\underset{\mathbf{p}}{\max }\mathrm{Re}⟨\mathbf{x},\mathbf{p}⟩\text{mbox}subjectto{\Vert \mathbf{p}\Vert }_{\mathcal{A}}^{\ast }\le 1,$$

$${\Vert \mathbf{p}\Vert }_{\mathcal{A}}^{\ast }\triangleq$$

$$\underset{\mathbf{p}\in {\mathbb{C}}^n,\mathbf{H}\in {\mathbb{C}}^{n\times n}}{\max }$$

$$\left[\begin{array}{cc}\mathbf{H}& \mathbf{p}\\ {\mathbf{p}}^{\mathsf{H}}& 1\end{array}\right]⪰0$$

$$\sum _{i=1}^{n-k}{H}_{i,i+k}={\delta }_k,k=0,\dots ,n-1,$$

$$\hat{\mathcal{T}}=\left\{\tau :|\hat{P}(\tau )|=1\right\}.$$

$$\mathrm{Re}⟨\mathbf{x},\hat{\mathbf{p}}⟩$$

$$P_{\star }({\tau }_k)$$

$$\left|P_{\star }(\tau )\right|$$

$${\Delta }_{\mathbb{T}}\left(\mathcal{T}\right)\triangleq \underset{\begin{array}{c}\tau ,{\tau }^{\prime }\in \mathcal{T}\\ \tau \ne {\tau }^{\prime }\end{array}}{\inf }\underset{q\in \mathbb{Z}}{\min }\left|\tau -{\tau }^{\prime }+q\right|.$$

$${\Delta }_{\mathbb{T}}\left(\mathcal{T}\right)>\frac{2.52}{n-1}.$$

$$\underset{\mathbf{x}}{\min }\frac{1}{2}{\Vert \mathbf{x}-\mathbf{z}\Vert }_2^2+\lambda {\Vert \mathbf{x}\Vert }_{\mathcal{A}},$$

$$\frac{1}{n}{\Vert \hat{\mathbf{x}}-\mathbf{x}\Vert }_2^2=O\left(\sigma \sqrt{\frac{\log n}{n}}\sum _{k=1}^r|c_k|\right)$$

$$\frac{1}{n}{\Vert \hat{\mathbf{x}}-\mathbf{x}\Vert }_2^2=O\left({\sigma }^2\frac{r\log n}{n}\right).$$

Full text

Harnessing Sparsity over the Continuum:

Atomic Norm Minimization for Super Resolution

Yuejie Chi, Senior Member, IEEE, Maxime Ferreira Da Costa, Member, IEEE The authors are with the Department of Electrical and Computer Engineering, Carnegie Mellon University, Pittsburgh, PA, USA (emails: {yuejiechi, mferreira}@cmu.edu).

At the core of many sensing and imaging applications, the signal of interest can be modeled as a linear superposition of translated or modulated versions of some template (e.g. a point spread function, a Green’s function), and the fundamental problem is to estimate the translation or modulation parameters (e.g. delays, locations, Dopplers) from noisy measurements. This problem is of central importance to not only target localization in radar and sonar, channel estimation in wireless communications, direction-of-arrival estimation in array signal processing, but also modern imaging modalities such as super-resolution single-molecule fluorescence microscopy, nuclear magnetic resonance imaging, spike localization in neural recordings, among others.

Typically, the temporal or spatial resolution of the acquired signal is limited by that of the sensing or imaging devices, due to factors such as the numerical aperture of a microscope, the wavelength of the impinging electromagnetic or optical waves, or the sampling rate of an analog-to-digital converter. This resolution limit is well-known and often referred to as the Rayleigh limit (c.f. the box “What is the Rayleigh Limit?”). The performance of matched filtering, or periodogram, which creates a correlation map of the acquired signal against the range of parameters, is limited by the Rayleigh limit, regardless of the noise level.

On the other hand, the desired resolution of parameter estimation can be much higher – a challenge known as super resolution. There is a long history of pursuing super-resolution algorithms in the community of signal processing [1, 2, 3]. The oldest one probably dates back to de Prony’s root finding method in as early as 1795 [4], and variants of this method that are better suited for noisy data have also been proposed over time, see e.g. [5]. Subspace methods based on the computation of eigenvector or singular vector decompositions, such as MUSIC [6], ESPRIT [7] and matrix pencil [8], are another class of popular approaches since their inception in 1980s. Different forms of maximum likelihood estimators have also been studied extensively [9, 10]. Collectively, these algorithms have “super resolution” capabilities, namely, they can resolve the parameters of interest at a resolution below the Rayleigh limit when the noise level is sufficiently small.

While there already exists a plethora of traditional methods, convex optimization recently emerges as a compelling framework for performing super resolution, garnering significant attention from multiple communities spanning signal processing, applied mathematics, and optimization. Due to (the relative) tractability of convex analysis and convex optimization, the new framework offers several benefits. First, strong theoretical guarantees are rigorously established to back up its performance even in the presence of noise and corruptions. Second, it is versatile to include prior knowledge into the convex program to handle a wide range of measurement models that fall out of the reach of traditional methods. Third, leveraging the rapid progress in large-scale convex optimization, it opens up the possibility of applying efficient solvers tailored to real-world applications.

The goal of this paper is to offer a friendly exposition to atomic norm minimization [11] as a canonical convex approach for super resolution. The atomic norm is first proposed in [12] as a general framework for designing tight convex relaxations to promote “simple” signal decompositions, where one seeks to use a minimal number of “atoms” to represent a given signal from an atomic set composed of an ensemble of signal atoms. Celebrated convex relaxations such as the $\ell_{1}$ norm approach for cardinality minimization [13] and the nuclear norm approach for rank minimization [14] can be viewed as particular instances of atomic norms for appropriately defined atomic sets. Specializing the atomic set to a dictionary containing all translations of the template signal over the continuous-valued parameter space, estimating the underlying translation parameters is then equivalent to identifying a sparse decomposition in an infinite-dimensional dictionary. This key observation allows one to recast super resolution as solving an infinite-dimensional convex program [15] – a special form of atomic norm minimization considered in this paper. We first highlight its mathematical formulation through a pedagogical yet useful model of super resolution that amounts to line spectrum estimation, where this infinite-dimensional convex program can be equivalently reformulated as a semidefinite program, then demonstrate its versatility by discussing how it can be adapted to address measurement models that traditional methods may not be easily applicable. Finally, we illustrate its utility in super resolution image reconstruction for single-molecule fluorescence microscopy [16], where the infinite-dimensional convex program can be solved efficiently via tailored solvers.

Throughout this paper, we use boldface letters to represent matrices and vectors, e.g. ${\boldsymbol{a}}$ and ${\boldsymbol{A}}$. We use ${\boldsymbol{A}}^{\top}$, ${\boldsymbol{A}}^{\mathsf{H}}$, $\mathop{\rm Tr}({\boldsymbol{A}})$ to represent the transpose, Hermitian transpose, and trace of ${\boldsymbol{A}}$, respectively. The conjugate of a complex scalar $a$ is denoted as $a^{*}$. We use ${\boldsymbol{A}}\succeq 0$ to represent ${\boldsymbol{A}}$ is positive semidefinite. The matrix $\mathop{\rm toep}\left({\boldsymbol{u}}\right)$ denotes the Hermitian Toeplitz matrix whose first column is equal to ${\boldsymbol{u}}$, and $\mbox{diag}(\bm{g})$ denotes the diagonal matrix with diagonal entries given as $\bm{g}$. The inner product between two matrices ${\boldsymbol{X}}$ and ${\boldsymbol{P}}$ is defined as $\left\langle{\boldsymbol{X}},{\boldsymbol{P}}\right\rangle=\mathop{\rm Tr}({\boldsymbol{X}}^{\mathsf{H}}{\boldsymbol{P}})$. Additionally, the notation $f(n)=O\left(g(n)\right)$ means that there exists a constant $c>0$ such that $\left|f(n)\right|\leq c|g(n)|$.

I What is the Atomic Norm?

An everlasting idea in signal processing is decomposing a signal into a linear combination of judiciously chosen basis vectors, and seeking compact and interpretable signal representations that are useful for downstream processing. For example, decomposing time series into sinusoids, speeches and images into wavelets, total system responses into impulse responses, etc.

To fix ideas, consider the task of representing a signal ${\boldsymbol{x}}$ in a vector space using atoms from a collection of vectors in $\mathcal{A}=\{{\boldsymbol{a}}_{i}\}$ called an atomic set. The set ${\mathcal{A}}$ can contain either finite or infinite numbers of atoms. We wish to expand ${\boldsymbol{x}}$ using the atoms in a form of

[TABLE]

where $c_{i}>0$ specifies the coefficients of the decomposition. In many applications, the size of ${\mathcal{A}}$ can be much larger than the dimension of the signal, leading to an overcomplete representation, and there are an infinite number of possibilities to decompose ${\boldsymbol{x}}$. Which representation, then, shall we pick? Among the many plausible criteria, one meaningful approach is to pursue the Occam’s razor principle, and seek for a parsimonious decomposition of the signal ${\boldsymbol{x}}$ involving the smallest possible number of atoms in ${\mathcal{A}}$, i.e. the sparsest solution to (1). The corresponding representation is known as a sparse representation [17]. Many real-world signals admit sparse representations for appropriately chosen atomic sets. As a simple example, natural images are approximately sparse by picking $\mathcal{A}$ as a wavelet frame. Low-rank matrices, another class of signals that have enjoyed wide success in signal processing [18], are sparse with respect to an atomic set $\mathcal{A}$ that is the collection of all unit-norm rank-one matrices.

Given a signal ${\boldsymbol{x}}$, how to find its sparse representation in the atomic set ${\mathcal{A}}$? In general, this problem is nonconvex and can be NP-hard due to the combinatorial aspect of cardinality minimization. The key motivation behind atomic norm minimization, proposed by Chandrasekaran et al. [12], is to relax the nonconvex sparsity cost by its tight convex surrogate, and solve instead the resulting convex relaxation that is more tractable. This idea is a generalization of the popular $\ell_{1}$ minimization for sparse vector recovery [19, 20] when ${\mathcal{A}}$ is a finite set. Therein, one seeks to solve a linear program which minimizes the sum instead of the cardinality of the nonzero coefficients.

To extend the same idea to the case where ${\mathcal{A}}$ is an arbitrary, and possibly infinite-dimensional set, we first take the convex hull of ${\mathcal{A}}$, denoted as $\mathop{\rm conv}({\mathcal{A}})$, and then define its associated Minkowski functional (or gauge function) as [12]

[TABLE]

which is the solution to a convex program. When ${\mathcal{A}}$ is centrally symmetric about the origin, the above definition leads to a valid norm, and is called the atomic norm of ${\boldsymbol{x}}$. Fig. 2 presents an illustration of this concept, where the atomic norm is the smallest nonnegative scaling of $\mathop{\rm conv}({\mathcal{A}})$ until it intersects ${\boldsymbol{x}}$. Following the definition (2), a fundamental geometric property is that the atomic norm ball, i.e., $\{{\boldsymbol{x}}:\;\|{\boldsymbol{x}}\|_{{\mathcal{A}}}\;\leq 1\}$, is exactly $\mathop{\rm conv}({\mathcal{A}})$.

More interestingly, consider the case when ${\boldsymbol{x}}$ lies in an $n$-dimensional vector space. Carathéodory’s theorem [21] guarantees that any point in $\mathop{\rm conv}({\mathcal{A}})$ can be decomposed as a convex combination of at most $n+1$ points in ${\mathcal{A}}$, where ${\mathcal{A}}$ is not necessarily convex. Therefore, one may rewrite (2) as

[TABLE]

as long as the centroid of $\mathop{\rm conv}({\mathcal{A}})$ is the origin. The decomposition $\sum_{i}c_{i}{\boldsymbol{a}}_{i}$ that obtains the infimum is referred to as the atomic decomposition of ${\boldsymbol{x}}$ onto ${\mathcal{A}}$. It is not hard to see that the atomic norm indeed subsumes the $\ell_{1}$ norm as a special case but accommodates the more general case where ${\mathcal{A}}$ can be an infinite-dimensional set.

Several central questions are how to properly select the atomic set, compute the atomic norm and find the atomic decomposition, and when the atomic decomposition coincides with the sparse representation, i.e. the convex relaxation is tight. Clearly, the answers depend highly on the atomic set as well as the signal itself. These questions have been addressed extensively in the study of $\ell_{1}$ norm minimization for sparse vector recovery [19, 20, 22]. In the context of super resolution, we will first address these questions under a simple model that amounts to the classical problem of line spectrum estimation, which has deep connections to systems and control theory.

II A Mathematical Model of Super Resolution, Equivalent to Line Spectrum Estimation

We first focus on a simple yet widely applicable model of super resolution that describes the convolution of a sequence of point sources with a point spread function (PSF) that is resolution-limited, illustrated in Fig. 3. Let $x(t)$ be a spike signal given as

[TABLE]

Here, $r$ is the number of spikes, $c_{k}\in\mathbb{C}$ and $\tau_{k}\in[0,1)$111Without loss of generality, the maximal delay is normalized to $1$. are the complex amplitude and delay of the $k$th spike. Such a spike signal can model many physical phenomena, such as firing times of neurons, locations of fluorescence molecules, and so on. Let $g(t)$ be the PSF whose bandwidth is limited due to the Rayleigh limit, namely its Fourier transform $G(f)$ satisfies

[TABLE]

for some bandwidth $B>0$. Its convolution with $x(t)$, contaminated by an additive noise $\epsilon(t)$, can be written as

[TABLE]

where $*$ denotes the convolution operator. Sampling the Fourier transform of the above equation at the frequencies $\ell=-\lfloor B/2\rfloor,\cdots,0,\cdots,\lfloor B/2\rfloor$, we obtain the measurements

[TABLE]

where $G_{\ell}$, $X_{\ell}$, $E_{\ell}$, and $Y_{\ell}$ are the Fourier transforms of $g(t)$, $x(t)$, $\epsilon(t)$, and $y(t)$ evaluated at frequency $\ell$, respectively. The total number of samples is $n=2\lfloor B/2\rfloor+1\approx B$. We write (5) in a vector form as

[TABLE]

where $\bm{y}=[Y_{\ell}]$, $\bm{g}=[G_{\ell}]$, $\bm{x}=[X_{\ell}]$, and $\bm{\epsilon}=[E_{\ell}]$. The problem of super resolution is then to estimate $\{c_{k},\tau_{k}\}_{1\leq k\leq r}$ accurately from $\bm{y}$, without knowing the model order $r$ a priori. Here, the Rayleigh limit is inversely proportional to the bandwidth $B$, and roughly speaking, is about $1/n$.

When the PSF $g(t)$ is known, one can “equalize” (5) by multiplying $G_{\ell}^{-1}$ to both sides, provided that $G_{\ell}$’s are non-zero. The observation ${\boldsymbol{z}}=[G_{\ell}^{-1}Y_{\ell}]$ relates to ${\boldsymbol{x}}$ as

[TABLE]

where $\widetilde{\bm{\epsilon}}$ is the additive noise. With slight abuse of notation, we map the index of $\ell$ from $\lfloor B/2\rfloor,\cdots,-\lfloor B/2\rfloor$ to $0,\cdots,n-1$ for convenience, and write $\bm{x}$ as a superposition of complex sinusoids:

[TABLE]

where ${\boldsymbol{a}}\left(\tau\right)\in\mathbb{C}^{n}$ is a vector defined as

[TABLE]

Notably, the above simplified model (7) also amounts to the classical problem of line spectrum estimation, that consists of estimating a mixture of sinusoids (with frequencies $\tau_{k}\in[0,1)$) from equi-spaced time samples (sampled at integers $\{0,\cdots,n-1\}$) of the time-domain signal $x_{\mathsf{ls}}(t)=\sum_{k=1}^{r}c_{k}e^{j2\pi\tau_{k}t}$. This finds applications in speech processing, power system monitoring, systems identification, and so on. The same model also describes direction-of-arrivals estimation using a uniform linear array, which is studied extensively in the literature of spectrum analysis [1].

III Line Spectrum Super Resolution via Atomic Norm Minimization

In the absence of noise, one could think of super resolution as estimating the continuous-time spike signal $x(t)$ in (4) from its discrete-time moment measurements ${\boldsymbol{x}}$ in (8), which are related through

[TABLE]

One can also think of $x(t)$ as the representation of $\bm{x}$ over a continuous dictionary

[TABLE]

which forms a one-dimensional variety of $\mathbb{C}^{n}$ called the moment curve, illustrated in Fig. 4 (a). It is well-known that the convex hull of ${\mathcal{A}}_{0}$, illustrated in Fig. 4 (b), is a body of $\mathbb{C}^{n}$ that can be parameterized by a set of linear matrix inequalities [23], and has close relationships with the positivity of Hermitian Toeplitz matrices. This fundamental property of the moment curve has many implications in control and signal processing [24, 25], and is key to the development of a super-resolution theory based on atomic norm minimization.

It is clearly possible to obtain the same ${\boldsymbol{x}}$ from different $x(t)$. However, if we impose some sparsity assumption, namely constraining how many spikes are allowed in $x(t)$, this representation can be ensured to be unique. In particular, the representation (8) is unique as long as $r\leq\left\lfloor n/2\right\rfloor$ and the support set $\mathcal{T}=\{\tau_{k}\}_{1\leq k\leq r}$ contains distinct elements.

III-A Atomic Norm for Line Spectrum Super Resolution

To apply the framework of atomic norm minimization for super resolution, one must first define the atomic set properly. Since the complex amplitudes $c_{k}$’s can take arbitrary phases, we introduce an augmented atomic set taking this into account:

[TABLE]

See an illustration of ${\mathcal{A}}_{\mathsf{1D}}$ and its convex hull in Fig. 4 (c) and (d). Writing $c_{k}=\left|c_{k}\right|e^{j\phi_{k}}$, ${\boldsymbol{x}}$ can be represented as a positive combination of the atoms in ${\mathcal{A}}_{\mathsf{1D}}$ as ${\boldsymbol{x}}=\sum_{k=1}^{r}\left|c_{k}\right|e^{j\phi_{k}}{\boldsymbol{a}}\left(\tau_{k}\right)$. It is easy to verify that ${\mathcal{A}}_{\mathsf{1D}}$ is centrally symmetric around the origin, and consequently it induces an atomic norm over $\mathbb{C}^{n}$, as defined in (2) and (3). It is worth noting that minimizing the atomic norm of $\bm{x}$ is equivalent to minimizing the total variation of $x(t)$, i.e.

[TABLE]

and both viewpoints are used frequently in the literature.

Remarkably, this atomic norm admits an equivalent semidefinite program (SDP) characterization, thanks to the Carathéodory–Fejér–Pisarenko decomposition [26]:

[TABLE]

Contrary to its abstract definition in (2), the reformulation (III-A) provides a tractable approach to compute the quantity $\left\|{\boldsymbol{x}}\right\|_{{\mathcal{A}}}$, which can be accomplished using generic off-the-shelf convex solvers [27]. The Vandermonde decomposition of $\mathop{\rm toep}({\boldsymbol{u}})$, i.e. $\mathop{\rm toep}({\boldsymbol{u}})=\sum_{k=1}^{r^{\prime}}|c_{k}^{\prime}|{\boldsymbol{a}}(\tau_{k}^{\prime}){\boldsymbol{a}}(\tau_{k}^{\prime})^{\sf{H}}$ can then be used to identify the support $\widehat{\mathcal{T}}=\{\tau_{k}^{\prime}\}$ of the atomic representation of ${\boldsymbol{x}}$, as well as the atomic norm $\|{\boldsymbol{x}}\|_{{\mathcal{A}}}=\sum_{k=1}^{r^{\prime}}|c_{k}^{\prime}|$.

III-B Duality and Atomic Decomposition

The Lagrangian duality theory marks an important aspect in understanding the atomic norm. The Lagrange dual problem associated with the atomic norm minimization (2) reads [11]

[TABLE]

where the dual atomic norm $\left\|{\boldsymbol{p}}\right\|_{{\mathcal{A}}}^{\ast}$ of a vector ${\boldsymbol{p}}\in\mathbb{C}^{n}$ is defined with respect to the atomic set $\mathcal{A}_{\mathsf{1D}}$ as

[TABLE]

The last equality of (17) suggests that the dual atomic norm can be interpreted as the supremum of the modulus of a complex trigonometric polynomial $P(\tau)=\left\langle{\boldsymbol{a}}\left(\tau\right),{\boldsymbol{p}}\right\rangle=\sum_{\ell=0}^{n-1}p_{\ell}e^{-j2\pi\ell\tau}$ with coefficients given by the vector ${\boldsymbol{p}}$. Constraints of this type is known to be equivalent to linear matrix inequalities involving the positivity of some Hermitian matrices (c.f. the box “From Bound Polynomials to Linear Matrix Inequalities”). The dual program (16) can be reformulated into the SDP below:

[TABLE]

where $H_{i,j}$ is the $(i,j)$th entry of the matrix ${\boldsymbol{H}}$, and the indicator function $\delta_{k}$ equals to $1$ if $k=0$ and [math] otherwise.

Another merit of the dual formulation is that the support set of the atomic decomposition can be inferred from the optimal solution $\widehat{{\boldsymbol{p}}}$ to the dual problem (16), by examining the dual polynomial $\widehat{P}(\tau)=\left\langle{\boldsymbol{a}}\left(\tau\right),\widehat{{\boldsymbol{p}}}\right\rangle$. We identify the spikes as the locations of the extreme values of the modulus of $\widehat{P}(\tau)$:

[TABLE]

This is possible, because, under strong duality, both the primal and the dual problems must share the same optimal objective value, i.e. $\left\|{\boldsymbol{x}}\right\|_{{\mathcal{A}}}=\sum_{k=1}^{r^{\prime}}|c_{k}^{\prime}|,$ where ${\boldsymbol{x}}=\sum_{k=1}^{r^{\prime}}c_{k}^{\prime}{\boldsymbol{a}}(\tau_{k}^{\prime})$ is the atomic decomposition of ${\boldsymbol{x}}$. Consequently, the optimal value of dual program becomes

[TABLE]

indicating $\widehat{P}(\tau_{k}^{\prime})=\mbox{sgn}(c_{k}^{\prime})=c_{k}^{\prime}/|c_{k}^{\prime}|$ whenever the atomic decomposition is non-vanishing at $\tau_{k}^{\prime}$. This approach is illustrated in Fig. 6 for a length-$33$ signal with 6 spikes, where the peaks of $\widehat{P}(\tau)$ matches with the locations of the true spikes, indicating the atomic decomposition perfectly recovers the true sparse representation.

The atomic norm offers an approach for line spectrum super resolution that is drastically different from traditional methods, which rely heavily on the correctness of model order estimation. The dual polynomial approach, in contrast, does not require any prior knowledge on the model order, and can estimate the spikes with an infinitesimal precision.

III-C Exact Recovery Guarantees

So far, we have explained the algorithmic approach of atomic norm minimization for line spectrum super resolution. A central question is to understand whether this convex relaxation is tight or not. More precisely, once would like to identify the conditions under which the estimated support $\widehat{\mathcal{T}}$ coincides with the true support $\mathcal{T}$ of the signal ${\boldsymbol{x}}$, and correspondingly, the atomic decomposition ${\boldsymbol{x}}=\sum_{k=1}^{r^{\prime}}c_{k}^{\prime}{\boldsymbol{a}}(\tau_{k}^{\prime})$ coincides with the sparsest representation ${\boldsymbol{x}}=\sum_{k=1}^{r}c_{k}{\boldsymbol{a}}(\tau_{k})$ over the atomic set ${\mathcal{A}}_{\mathsf{1D}}$.

Such questions were extensively addressed in the context of $\ell_{1}$ norm minimization, where the atomic set ${\mathcal{A}}$ has a finite number of elements. The performance guarantees often depend on specific structural properties of ${\mathcal{A}}$, formalized in the notion of restricted isometry property (RIP) [30], or certain incoherence properties [31]. Unfortunately, these properties do not hold when considering a continuous dictionary such as ${\mathcal{A}}_{\mathsf{1D}}$, since two atoms ${\boldsymbol{a}}\left(\tau\right)$ and ${\boldsymbol{a}}(\tau+\delta)$ can be more and more correlated with each other as their separation $\delta$ tends to zero, leading to arbitrarily small RIP or coherence constants.

Nonetheless, one could ask for which class of signals the relaxation is tight. Leveraging duality theory, the atomic norm approach is tight for a fixed signal ${\boldsymbol{x}}$, i.e. $\widehat{\mathcal{T}}=\mathcal{T}$, as long as there exists a dual certificate ${\boldsymbol{p}}_{\star}$, such that $P_{\star}(\tau)=\langle{\boldsymbol{a}}(\tau),{\boldsymbol{p}}_{\star}\rangle$ satisfies [11]

[TABLE]

In other words, it amounts to find an $(n-1)$-order trigonometric polynomial that interpolates sign patterns of the spike signal at the spike locations, as well as is bounded in magnitude by 1.

Intuitively, the difficulty of interpolation depends on the separations between the spikes in $\mathcal{T}$, and more precisely on the minimal separation, or the minimal wrap-around distance between any pair of distinct spikes in $\mathcal{T}$, defined formally as

[TABLE]

This metric is illustrated in Fig. 7, and reflects the periodic behavior of the atom ${\boldsymbol{a}}\left(\tau+q\right)={\boldsymbol{a}}\left(\tau\right)$ for any integer $q\in\mathbb{Z}$. For instance, if $\mathcal{T}=\{0.1,0.9\}$, then $\Delta_{\mathbb{T}}\left(\mathcal{T}\right)=0.2$.

A remarkable result, established by Candès and Fernandez-Granda in [11], is that for sufficiently large $n$, a valid certificate can be constructed in a deterministic fashion, as long as the separation condition $\Delta_{\mathbb{T}}\left(\mathcal{T}\right)\geq\frac{4}{n-1}$ holds, regardless of the complex amplitudes of the spikes. Furthermore, this result does not make any randomness assumptions on the signal. Later, this separation condition has been further improved by Fernandez-Granda [32] to

[TABLE]

Conversely, there exist some spike signals with $\Delta_{\mathbb{T}}\left(\mathcal{T}\right)<\frac{2}{n-1}$ such that atomic norm minimization fails to resolve [33].

III-D Atomic Norm Denoising

In practice, the observations are corrupted by noise, and no estimator can exactly recover the spike signal $x(t)$. This raises a natural question regarding the robustness of the estimate produced by atomic norm minimization methods. When the noise is additive and the observation ${\boldsymbol{z}}$ obeys the noisy model (7), it has been proposed to estimate ${\boldsymbol{x}}$ by searching around the observation ${\boldsymbol{z}}$ for signals with small atomic norms [35]:

[TABLE]

where $\lambda>0$ is a regularization parameter that draws a trade-off between the fidelity to the observation and the size of the atomic norm. This method, known as “atomic norm denoising”, can be interpreted as a generalization of the celebrated LASSO estimator [36].

When the noise vector $\widetilde{\bm{\epsilon}}$ is composed of i.i.d. complex Gaussian entries with zero mean and variance $\sigma^{2}$, the mean squared error (MSE) of the estimate $\widehat{{\boldsymbol{x}}}$ returned by (22) can be bounded as [35]

[TABLE]

with high probability by setting $\lambda=\eta\sigma\sqrt{n\log n}$ for some constant $\eta\in(1,\infty)$, e.g. $\eta=1.2$ in practice. This error rate can be significantly improved when the spikes satisfy the separation condition $\Delta_{\mathbb{T}}\left(\mathcal{T}\right)\geq\frac{4}{n-1}$, where with high probability one has [37]

[TABLE]

This last error rate is near-optimal up to some logarithmic factor, since no estimator can achieve an MSE below the rate $O\left(\sigma^{2}\frac{r\log\left(n/r\right)}{n}\right)$ [37].

A more important performance criteria in super resolution concerns the stability of the support estimate $\widehat{\mathcal{T}}$, which has been studied in [37]–[41]. When the spikes satisfy the separation condition $\Delta_{\mathbb{T}}\left(\mathcal{T}\right)>\frac{5.0018}{n-1}$, and the complex amplitudes of the coefficients $\left\{c_{k}\right\}_{1\leq k\leq r}$ have approximately the same modulus, then it is established in [39] that the atomic decomposition of the output $\widehat{{\boldsymbol{x}}}$ of (22) is composed of the same number of spikes, i.e. $|\widehat{\mathcal{T}}|=|\mathcal{T}|=r$ and that the estimated parameters satisfy

[TABLE]

with high probability. Altogether, it can be seen that atomic norm denoising achieves near-optimal performance guarantees as long as the spikes are separated by a few times the Rayleigh limit (c.f. the box “Is the Separation Condition Necessary?”).

It is natural to wonder how atomic norm denoising fares compared with classical approaches such as Prony and MUSIC for line spectrum estimation. We examine their ability to resolve close-located spikes with opposite signs, where the reconstruction performance is measured in terms of the MSE of the estimated spike locations $\widehat{\mathcal{T}}$, and for different values of separation $\Delta_{\mathbb{T}}(\mathcal{T})=\alpha/n$. Fig. 9 shows the MSE of atomic norm denoising, Prony’s method with Cadzow denoising [42], and root-MUSIC [43] with respect to the SNR defined as $\|{\boldsymbol{x}}\|_{2}^{2}/(n\sigma^{2})$, benchmarked against the Cramér-Rao bound (CRB), when the separation parameter $\alpha=2$, $1.75$, $1.5$ respectively. It is clear that atomic norm denoising outperforms classical approaches, and approaches the CRB at a much lower SNR.

III-E A Faster Algorithm via ADMM

While the SDP formulation is tractable, its computational complexity is prohibitive when solving large-dimensional problems. Fortunately, it is possible to develop tailored algorithms that are significantly faster. For conciseness, we will discuss one approach based on the Alternating Direction Method of Multipliers (ADMM) [35]. The general principle of ADMM is to split the quadratically-augmented Lagrangian function of an optimization problem into a sum of separable sub-functions [44]. Each iteration of the algorithm consists of performing independent local minimization on each of those quantities, while ensuring that the feasibility constraints are always satisfied. The iterations run until both primal and dual residuals satisfy a pre-defined tolerance level.

We take atomic norm denoising (22) as an example, which, in light of (III-A), can be equivalently rewritten as

[TABLE]

The above program has been “augmented” by introducing an intermediate variable $\bm{S}$ for the purpose of decoupling the positive semidefinite constraint on the matrix $\bm{S}$ from the linear constraints on its structure. The augmented Lagrangian $\mathcal{L}$ is given as

[TABLE]

where $\boldsymbol{S}$ and $\boldsymbol{\Sigma}$ are $(n+1)$-dimensional Hermitian matrices, and $\rho>0$ is a regularization parameter. The successive update steps to minimize the augmented Lagrangian are given in Alg. 1. Closed-form solutions can be found for the first update step, yielding a very efficient implementation. The second update is the most costly part, as a projection over the cone of positive semidefinite Hermitian matrices has to be computed. This computation is typically achieved using power methods [45], with a computational complexity of $O(n^{3})$ per iteration.

III-F Can we discretize?

It may be worthwhile to pause and compare atomic norm minimization to other approaches based on convex optimization for super resolution, in particular, $\ell_{1}$ minimization that is widely popular for high-resolution imaging and localization in the recent literature due to Compressed Sensing (CS) [46, 47].

The $\ell_{1}$ norm can be seen as a discrete approximation of the atomic norm. Indeed, taking the atomic set ${\mathcal{A}}_{\mathsf{1D}}$, one can pick a desired resolution $Q$ and discretize it as:

[TABLE]

and then perform $\ell_{1}$ minimization over ${\mathcal{A}}_{\mathsf{1D,discrete}}$. The convex hull of ${\mathcal{A}}_{\mathsf{1D,discrete}}$ approaches that of ${\mathcal{A}}_{\mathsf{1D}}$ as the discretization gets finer, which suggests the performance of $\ell_{1}$ minimization over the discretized dictionary approaches that of atomic norm minimization asymptotically [48]. If the spike signal meets a so-called “non-degenerate source condition” [49, Definition 2], this approach will return a sparse solution supported on the elements of the discretized ${\mathcal{A}}_{\mathsf{1D,discrete}}$ surrounding the ground truth spikes, when the noise is small enough [49, 50].222However, it remains unclear which class of spike signals satisfies the non-degenerate source condition in practice.

However, this discretization may come with several undesired consequences when the grid size $Q$ is finite in practice. The theory of $\ell_{1}$ minimization only provides exact recovery guarantees when the spikes of $x(t)$ lie on the grid, which is unrealistic. In fact, there is always an inevitable basis mismatch [51], between the spikes represented in the discretized dictionary ${\mathcal{A}}_{\mathsf{1D,discrete}}$ and the true spikes, no matter how fine the grid is. Perfect recovery is not possible in this situation even in the absence of noise due to this mismatch. Furthermore, one can find signals whose representations in ${\mathcal{A}}_{\mathsf{1D,discrete}}$ are not compressible due to spectral leakage, and therefore are poorly recovered using $\ell_{1}$ minimization, e.g. the recovery may contain many spurious spikes. Therefore, cautions are needed to account for such consequences when applying discretization, and efforts to mitigate the basis mismatch have been proposed extensively, e.g. [52, 53].

IV Generalizations of Atomic Sets

The tool of atomic norms can be extended easily to handle a wide range of scenarios in a unified manner, by properly adjusting the atomic set for signal decompositions, such as incorporating prior information, dealing with multi-dimensional settings and multiple measurement vectors, to illustrate a few.

IV-A Atomic Set for Positive Spikes

In some applications, there exist additional information about the spikes, such as the coefficients of the spikes in (8) are positive, i.e. $c_{k}>0$. Examples include neural spike sorting, fluorescence microscopy imaging, or covariance-based spectrum estimation for noncoherent sources [1].

In this case, the atomic set reduces to the moment curve ${\mathcal{A}}_{0}$ in (11). The induced $\|{\boldsymbol{x}}\|_{{\mathcal{A}}}$ is no longer a norm, since ${\mathcal{A}}_{0}$ is not centrally symmetric, but nonetheless, similar SDP characterization still holds. To be specific, the dual program now becomes

[TABLE]

where the constraint bounds the real part of the trigonometric polynomial $P(\tau)=\left\langle{\boldsymbol{a}}\left(\tau\right),{\boldsymbol{p}}\right\rangle$. Using the Fejér-Riesz Theorem (see e.g. [25, Theorem 1.1] and the box “From Bounded Polynomials to Linear Matrix Inequalities”), this can be equivalently represented as:

[TABLE]

It is long established [54, 55] that the spikes can be perfectly localized as long as $r\leq\lfloor(n-1)/2\rfloor$, without requiring any separation between the spikes, as long as they are positive. The stability of this approach as well as the implications of nonnegative constraints for other atomic sets are further studied in [56, 57, 58, 59]. As a comparison, Fig. 10 shows the MSE of atomic norm denoising with and without positive constraints, Prony’s method with Cadzow denoising [42], and root-MUSIC [43] with respect to the SNR defined as $\|{\boldsymbol{x}}\|_{2}^{2}/(n\sigma^{2})$ for resolving two spikes with positive signs, separated by $\Delta_{\mathbb{T}}(\mathcal{T})=\alpha/n$ for $\alpha=1,0.75,0.5$ respectively. It can be seen that atomic norm denoising still outperforms classical approaches, and in particular, incorporating the positive constraint leads to further improvements.

IV-B Atomic Set for Multi-Dimensional Spikes

When the spikes reside in a multi-dimensional space, one can extend the one-dimensional model in a straightforward manner. Here, we illustrate the setup for the two-dimensional case, where each entry of the signal $\boldsymbol{X}_{\mathsf{2D}}\in\mathbb{C}^{n_{1}\times n_{2}}$ can be expressed as a superposition of $r$ complex sinusoids propagating in two directions:

[TABLE]

where $c_{k}$ and $\bm{\tau}_{k}=\left[\tau_{1,k},\tau_{2,k}\right]^{\mathsf{T}}\in[0,1)^{2}$ are the complex amplitude and location of the $k$th spike, and $\boldsymbol{a}_{i}(\tau)$ is given by (9) with the dimension parameter replaced by $n_{i}$, $i=1,2$. It is natural to define the corresponding atomic set as [11, 60]

[TABLE]

and the atomic norm according to (2). To localize the spikes, one could similarly study the associated dual problem,

[TABLE]

where the dual atomic norm can be reinterpreted as the supremum of a bivariate complex trigonometric polynomial $P\left(\tau_{1},\tau_{2}\right)=\langle{\boldsymbol{a}}_{1}(\tau_{1}){\boldsymbol{a}}_{2}(\tau_{2})^{\top},{\boldsymbol{P}}\rangle$ with the matrix ${\boldsymbol{P}}$ as its coefficients. Again, one can localize the spikes by examining the extremal points of the dual polynomial, which is illustrated in Fig. 11. Cautions need to be taken when attempting to solve the dual program in two or higher dimensions, since the Bounded Real Lemma [24, 25] does not hold anymore. Instead, a precise characterization requires solving a hierarchy of sum-of-squares relaxations, and fortunately in practice, the first level usually suffices [61, 24, 60].

The tightness of the atomic norm minimization approach is closely related to a separation condition analogous to the one-dimensional case [32]. Namely, the atomic decomposition is unique and exact, as soon as there exists a universal constant $C>0$ such that the set of spikes $\mathcal{T}=\{\bm{\tau}_{k}\}_{1\leq k\leq r}$ satisfies

[TABLE]

Moreover, if the signal $\bm{X}_{\mathsf{2D}}$ is real-valued, $C=4.76$ suffices to guarantee exact recovery of the spikes.

IV-C Atomic Set for Multiple Measurement Vectors

One can collect multiple snapshots of observations, where they share the same spike locations with varying coefficients. Consider $T$ snapshots, stacked in a matrix, ${\boldsymbol{X}}_{\mathsf{MMV}}=[{\boldsymbol{x}}_{1},\ldots,{\boldsymbol{x}}_{T}]$, which is expressed similarly to (8) as

[TABLE]

where ${\boldsymbol{c}}_{k}=[c_{1,k},\ldots,c_{L,k}]\in\mathbb{C}^{T}$ is the coefficient of the $k$th spike across the snapshots. Following the recipe of atomic norms, we define the atoms as

[TABLE]

where $\tau\in[0,1)$, ${\boldsymbol{b}}\in\mathbb{C}^{T}$ with $\|{\boldsymbol{b}}\|_{2}=1$. The atomic set is defined as

[TABLE]

The atomic norm can be then defined following (2) which turns out sharing similar nice SDP characterizations for primal and dual formulations as for the single snapshot model [62]. The atomic norm $\|{\boldsymbol{X}}_{\mathsf{MMV}}\|_{\mathcal{A}}$ can be written equivalently as

[TABLE]

A curious comparison can be drawn to the nuclear norm by noticing that one recovers the nuclear norm of ${\boldsymbol{X}}_{\mathsf{MMV}}$ by replacing the principal block $\mathop{\rm toep}({\boldsymbol{u}})$ in (IV-C) with an arbitrary positive semidefinite matrix. The fact that $\mathop{\rm toep}({\boldsymbol{u}})$ has significantly fewer degrees of freedom ($n$ versus $n^{2}$) is in parallel to that ${\boldsymbol{a}}(\tau)$ has significantly fewer degrees of freedom than an arbitrary vector ($1$ versus $n$).

Again, one can determine the atomic decomposition and localize the spikes by resorting to the dual program in a similar fashion. Fig. 12 illustrates an example that multiple snapshots improve the performance of localization over the single snapshot case when the coefficients across snapshots exhibit some kind of diversity, e.g. generated with i.i.d. complex Gaussian entries.

V Generalizations of Measurement Models

So far, we have seen that atomic norm minimization provides a means for super resolution via convex relaxation in additive Gaussian noise. The framework of convex optimization is quite versatile and can be extended to handle models when the measurements are partially observed, corrupted, contain interfering sources, or even come from unknown modulations. This is an important advantage over classical methods such as MUSIC or ESPRIT, as most of them cannot be extended easily to these variants of models.

V-A Compressed spectral sensing

CS [46, 47] has suggested that it is possible to recover a signal using a number of measurements that is proportional to its degrees of freedom, rather than its ambient dimension. Consider the problem where only a subset of entries of $\bm{x}$ is observed,

[TABLE]

where $\bm{A}_{\mathsf{CS}}\in\mathbb{C}^{m\times n}$, and $m\ll n$ representing compressive acquisition of the signal ${\boldsymbol{x}}$. The goal is to recover $\bm{x}$ and its spectral content from $\bm{y}_{\mathsf{CS}}\in\mathbb{C}^{m}$, the compressive measurements. This has applications in wideband spectrum sensing and cognitive radio [63], for example.

One can easily extend the framework of atomic norm minimization, and recover $\bm{x}$ by solving the following program

[TABLE]

When $\bm{A}_{\mathsf{CS}}$ is a partial observation matrix, namely, a subset of $m$ entries of ${\boldsymbol{x}}$ is observed uniformly at random, then ${\boldsymbol{x}}$ can be perfectly recovered with high probability using $m=O\left(\log^{2}n+r\log r\log n\right)$ measurements as long as ${\boldsymbol{x}}$ satisfies the separation condition, and with random signs of the coefficients $c_{k}$’s [15, Theorem II.3]. More generally, a broader class of measurement matrices ${\boldsymbol{A}}_{\mathsf{CS}}$ can be allowed where its rows are drawn independently from some isotropic and incoherent distribution [64, 65], and exact recovery is possible under the same separation condition using a number of measurements on the order of $r$ up to some logarithmic factors. In addition, quantized measurements are further dealt in [66] with theoretical guarantees.

V-B Demixing sinusoids and spikes

Due to sensor failures or malicious environments, the measurements are susceptible to corruptions that can take arbitrary magnitudes. To this end, consider the problem when the observations are contaminated by sparse outliers, where

[TABLE]

Here, $\bm{s}$ is a sparse vector, where its nonzero entries correspond to corruptions of the observations. The goal is to decompose $\bm{x}$ and $\bm{s}$ from $\bm{y}_{\mathsf{corrupt}}$, a problem intimately related to the uncertainty principle of signal decomposition in [67, 17] and sparse error correction in CS [68].

Leveraging low-dimensional structures in both ${\boldsymbol{x}}$ and ${\boldsymbol{s}}$, we seek $\bm{x}$ with a small atomic norm and $\bm{s}$ with a small $\ell_{1}$ norm that satisfies the observation constraint [69]:

[TABLE]

where $\lambda>0$ is some regularization parameter. As long as the sample size is sufficiently large [69, Theorem 2.2], and the spikes satisfy the separation condition, then the above algorithm perfectly localizes the spikes with high probability, even when the corruption amounts to near a constant fraction of the measurements.

V-C Demixing interfering sources

A scenario of increasing interest is when the observation is composed of a mixture of responses from multiple exciting or transmitting sources, and the goal is to simultaneously separate and localize the sources at a high resolution. For example, an electrode probing the activities in a brain records firing patterns of a few neighboring neurons, each with a distinct PSF. For pedagogical reasons, let us consider a generalization of the model (6) with two interfering sources, where the observation is given as

[TABLE]

where $\bm{g}_{1}$ and $\bm{g}_{2}$ correspond to the frequency-domain response of the PSFs, and $\bm{x}_{i}=\sum_{k=1}^{r_{i}}c_{i,k}\bm{a}(\tau_{i,k})$, $i=1,2$. The goal is to separate and recover the spikes in both $\bm{x}_{1}$ and $\bm{x}_{2}$ from $\bm{y}_{\mathsf{mix}}$, where $\bm{g}_{1}$ and $\bm{g}_{2}$ are assumed known.

Using atomic norm minimization, one seeks to recover both ${\boldsymbol{x}}_{1}$ and ${\boldsymbol{x}}_{2}$ simultaneously by minimizing the weighted sum of their atomic norms [70]:

[TABLE]

where $\lambda>0$ is some regularization parameter. Unlike the single source case, the success of demixing critically depends on how easy it is to tell two PSFs apart – the more similar $\bm{g}_{1}$ and $\bm{g}_{2}$ are, the harder it is to separate them. A random model can be used to generate dissimilar PSFs, namely, it is assumed the entries of $g_{i}$’s are i.i.d. generated from a uniform distribution over the complex circle [70]. The algorithm then succeeds with high probability as long as the sample size is sufficiently large and the spikes within the same signal satisfy the separation condition [70, Theorem 2.1], without requiring any separation for spikes coming from different sources.

V-D Blind super resolution

So far, all algorithms have assumed the PSF as known, which is a reasonable assumption for problems where one can design and calibrate the PSF a priori. In general, one might need to estimate the PSF at the same time, possibly due to the fact that the PSF may drift and needs to be calibrated on the fly during deployment. In this case, we need to revisit (6) and estimate simultaneously $\bm{g}$ and $\bm{x}$ from their bilinear measurements,

[TABLE]

This problem is terribly ill-posed, as the number of unknowns far exceeds the number of observations. One remedy is to exploit additional structures of $\bm{g}$. For example, if $\bm{g}$ lies in a known low-dimensional subspace $\bm{B}=[\bm{b}_{1},\cdots,\bm{b}_{n}]^{\top}\in\mathbb{C}^{n\times d}$ with $d\ll n$, then the degrees of freedom of ${\boldsymbol{g}}$ is greatly dropped, since one only needs to estimate the coefficient $\bm{h}\in\mathbb{C}^{d}$ of ${\boldsymbol{g}}=\bm{B}\bm{h}$ in that subspace, which is of much smaller dimension. Even such, the measurements $\bm{y}_{\mathsf{BR}}$ is still bilinear in $\bm{h}$ and $\bm{x}$, and one cannot directly apply atomic norm minimization to ${\boldsymbol{x}}$ as it does not lead to a convex program.

Interestingly, a lifting trick can be applied [71], which rewrites $\bm{y}_{\mathsf{BR}}=\mathcal{X}(\bm{Z})$ as linear measurements of a higher-dimensional object $\bm{Z}=\bm{x}\bm{h}^{\top}\in\mathbb{C}^{n\times d}$ similar to (24):

[TABLE]

where $\bm{e}_{i}$ is the $i$th standard basis vector. Consequently, one can apply atomic norm minimization to $\bm{Z}$ with respect to (IV-C), leading to the algorithm [71]:

[TABLE]

This approach succeeds with high probability as soon as the sample size is sufficiently large, the spikes are well-separated and the PSF satisfies certain incoherence properties [71, Theorem 1]. Moreover, it can be further extended to demixing a mixture of sources with unknown PSFs, where each PSF lies in a distinct subspace, see [72].

VI Beyond Line Spectrum Estimation: Super-Resolution Imaging for Single-Molecule Fluorescence Microscopy

When the atomic set is composed of a family of complex sinusoidal signals, an exact implementation of the atomic norm in SDP is available. In the most general setting, atomic norm minimization is an infinite-dimensional convex program whose computation needs to be addressed carefully. Encouragingly, tailored solvers have been proposed and applied successfully to practical applications such as super-resolution imaging for single-molecule fluorescence microscopy, which we now present as a case study to show its promise.

VI-A Imaging Principle

The development of super-resolution fluorescence microscopy, which has been awarded the 2014 Nobel prize in Chemistry, is considered to fundamentally impact biological science and medicine. To date, a partial list of super-resolution fluorescence microscopy technologies includes PALM [73], STORM [74], and fPALM [75] which share a similar imaging principle. A very nice introduction can be found at [76]. While optical microscopy is desirable for imaging complex biological processes in live cells due to its noninvasive nature, due to diffraction limit, which is about hundreds of nanometers, it cannot image detailed internal structures of cells, which are often below 100 nanometers.

To deal with this challenge, biologists have come up with a clever idea of divide-and-conquer. To begin, imagine that every point within a cell is equipped with a photoswitchable fluorescent molecule, which means, once excited, the molecule will emit light stochastically over a duration of time to identify its location. This allows one to divide the imaging process into many frames, where in each frame, a random and sparse subset of fluorescent molecules (point sources) are activated and localized at a resolution below the diffraction limit using imaging algorithms. The final image is thus obtained by superimposing the localization outcomes of all the frames. Therefore, the high spatial resolution is achieved by sacrificing the temporal resolution. To speed up the imaging process and improve the temporal resolution, it is desirable to develop localization algorithms that are capable of identifying more fluorescent molecules per frame, which is known as the emitter density.

Very interestingly, this imaging principle can be used to reconstruct a 3-D biological structure from 2-D image frames [77]. One way is to introduce a cylindrical lens to modulate the ellipticity of the PSF based on the depth of the fluorescent object, which can be modeled as a Gaussian pulse with varying ellipticity along the $x$ and $y$ directions,

[TABLE]

where $\sigma_{x}(z)$ and $\sigma_{y}(z)$ are functions of the depth in the $z$ direction, and can be calibrated in advance. For a 3-D scene of point sources,

[TABLE]

its “convolution” with the PSF $g(x,y|z)$ is given as a 2-D image in the form of

[TABLE]

Therefore, to perform super-resolution, one needs to decode simultaneously the ellipticity as well as the location of the PSF, which is much more challenging. In practice, the situation is even more complex, since the continuous spatial function $(\zeta*g)(x,y)$ needs to be further discretized due to pixelization of the detector, leading to a discretized 2-D image, ${\boldsymbol{w}}\in\mathbb{R}^{n_{1}\times n_{2}}$, where each entry of ${\boldsymbol{w}}$ corresponds to the integration of $(\zeta*g)(x,y)$ over the area of a pixel. The final image, ${\boldsymbol{z}}$, counting the number of photons hitting the detector at each pixel, is modeled as a Poisson distribution with rate ${\boldsymbol{w}}$. The whole process is summarized in Fig. 13.

VI-B Applying Atomic Norm Minimization

Luckily, due to linearity, ${\boldsymbol{w}}$ can be viewed a sparse superposition of atoms that are parameterized by the 3-D point sources,

[TABLE]

where each atom ${\boldsymbol{a}}(x,y,z)$ corresponds to the image of a point source at $(x,y,z)$ after convolution and pixelization. The atomic set is then given as ${\mathcal{A}}_{\mathsf{3D}}=\{{\boldsymbol{a}}(x,y,z):x,y,z\in\mbox{imaging range}\}$. The goal is thus to recover $\zeta(x,y,z)$, or the atomic decomposition of ${\boldsymbol{w}}$, from the observation ${\boldsymbol{z}}\sim\mathsf{Poisson}({\boldsymbol{w}})$ as accurately as possible.

A natural approach is to seek the sparsest ${\boldsymbol{w}}$ such that the likelihood function of the observation ${\boldsymbol{z}}$ is maximized. To that end, we consider a constrained maximum likelihood estimation, where we seek to solve ${\boldsymbol{w}}$ via

[TABLE]

where $p({\boldsymbol{z}}|{\boldsymbol{w}})$ is the Poisson likelihood function, $\|{\boldsymbol{w}}\|_{{\mathcal{A}}}$ is the induced atomic norm with respect to ${\mathcal{A}}_{\mathsf{3D}}$, and $\eta$ is some regularization parameter that may be tuned in practice.

Early efforts such as CSSTORM [78], which are based on $\ell_{1}$ minimization by directly discretizing the parameter space, suffer from high computational complexity, due to the need of storing and manipulating a large dictionary of atoms, as fine discretization is required along all three spatial dimensions. On the other end, recent algorithmic developments such as ADCG [79] and CoGEnT [80] solve sparse inverse problems over continuous dictionaries with general convex loss functions at much reduced memory and computation requirements. In a nutshell, the ADCG method is an acceleration of conditional gradient, also known as Frank-Wolfe, to solve (26) with a general atomic set ${\mathcal{A}}=\{{\boldsymbol{a}}(\theta):\theta\in\Theta\}$, where $\theta$ is a short-hand notation for the parameter space. In particular, it directly estimates the atomic decomposition of ${\boldsymbol{w}}=\int{\boldsymbol{a}}(\theta)d\zeta(\theta)$.

The standard Frank-Wolfe adds one new atom at every iteration to reduce the negative log-likelihood function, however it will introduce many spurious atoms and lose sparsity as the iteration increases. To deal with this, ADCG introduces pruning and local refinements with a hope to maintain a sparse representation at all iterations. The detailed procedure of ADCG is given in Alg. 2. In the $j$th iteration, denote the support and coefficient of the current estimate of $\zeta(\theta)$ as $\mathcal{T}_{j}$ and ${\boldsymbol{c}}_{j}$, and the current estimate of ${\boldsymbol{w}}$ as $\mathcal{A}(\mathcal{T}_{j}){\boldsymbol{c}}_{j}$. Like Frank-Wolfe, ADCG starts by adding a spike to the estimated support $\mathcal{T}_{j}$ that maximally correlates with the derivative of $\log p({\boldsymbol{z}}|{\boldsymbol{w}})$ with respect to ${\boldsymbol{w}}$ at the current estimate, $\mathcal{A}(\mathcal{T}_{j})\bm{c}_{j}$. Because the spike location will be refined next, in practice, this step can be solved approximately by searching over a coarse grid of $\Theta$ to save computation.

ADCG then deviates from the standard Frank-Wolfe, and tries to improve the updated estimate by performing alternating descent over the coefficient and the support. It iterates between coefficients update via $\ell_{1}$ minimization, support pruning, and local refinement of the support by holding the coefficients fixed. The last step leverages the fact that $\bm{a}(\theta)$ is differentiable with respect to $\theta$, and a simple local search via gradient descent allows one to adjust the support to further reduce the loss function. Despite the nonconvexity in this refine step, [79] guarantees the convergence of ADCG with a convergence rate of $O(1/\epsilon)$ to reach $\epsilon$-accuracy in the function value, under some technical assumptions. In practice, the main computational benefits of ADCG are the absence of semidefinite constraints and small memory footprints, making it highly suitable for large-scale implementations.

TVSTORM [81] is a modification of ADCG tailored to solve (26) for 3-D image reconstruction with some domain adaptations to speed up implementations. TVSTORM outperforms CSSTORM both in terms of computational time and reconstruction quality. Fig. 14 (a) shows the diffraction-limited imaging using conventional microscopy; in contrast, the 3-D super-resolution image reconstructed using TVSTORM in Fig. 14 (b) is much clearer, where the structure of 3D microtubules can be well resolved with the axial coordinate represented in different colors. Fig. 14 (c) and (d) compare the reconstruction quality of a zoom-in region between TVSTORM and CSSTORM, where TVSTORM provides a visually more smooth reconstruction of the line structures in microtubules. Fig. 15 shows that TVSTORM indeed has a higher detection rate and a lower false discovery rate than CSSTORM, while executes much faster.

VII Final Remarks

In this paper, we presented an overview on how to leverage sparsity for continuous parameter estimation via the mathematical concept of atomic norms, which can be regarded as a generalization of the principle of $\ell_{1}$ norms for discrete model selection. We showcased its application in super resolution from low-pass observations in single-molecule fluorescence microscopy. The appeals of the atomic norm approach stem from its elegant mathematical framework, strong performance guarantees, and promises to scalable numerical solvers.

The atomic norm is only one of many possible approaches to exploit sparsity over the continuum. One competitive alternative is structured low-rank matrix optimization [82, 83, 84, 85]; see [86] for its connections and comparisons with the atomic norm approach. Another line of work [87, 88, 89] generalizes the traditional CS to an infinite-dimensional Hilbert space. In addition, sampling theorems are developed for signals with a finite rate of innovations together with strategies for perfect reconstruction [90]. More recently, a sparse functional framework has been proposed as a variational approach to handle sparsity over continuous and possibly nonlinear dictionaries. This category of estimators aims to recover functions with minimum support measure subject to the observation constraint [91, 92].

As a topic still under development, open problems abound for both theoretical and practical performance of optimization-based super resolution in general. We conclude by outlining some exciting future directions.

• •

Tight performance analysis in noise. Existing analyses of atomic norm denoising (22) typically only produce bounds that are tight up to some constant, making it less useful in practice. For example, attempts to benchmarking the theoretical bounds against the Cramér-Rao bound will be in vein due to the presence of large constants. It is therefore desirable to obtain tight performance bounds such as the one available for matrix denoising [93] that is asymptotically exact.

• •

Adaptive selection of regularization parameters. One benefit of atomic norm minimization over traditional spectrum estimation approaches is that it can automatically select the model order. However, the choice of the regularization parameter (22) depends on the noise level, which may not be available in practice. How to optimally set the regularization parameters is another problem of great importance, see [94] for some recent development.

• •

Low-rank factorization for SDP formulations of atomic norms. A popular heuristic to SDP with low-rank solutions is to apply low-rank matrix factorization and solve the corresponding nonconvex optimization problem [95, 96], with the premise of greatly reducing its computational cost. It will be interesting to see if this approach can be applied to speed up the computation of atomic norms with performance guarantees.

• •

Bridging classical and optimization-based approaches. There are deep connections between traditional approaches (e.g. Prony, MUSIC, and so on) and optimization-based approaches (e.g. atomic norm minimization and nuclear norm minimization). Such connections have already been realized, for instance, in the early works of Fuchs [97], where he provided an optimization interpretation of the Pisarenko method [98]. Another recent work [99] provided an optimization view to the MUSIC algorithm. It is hopeful that a confluence of past and current ideas will likely deepen our understandings and lead to further algorithmic improvements.

• •

Atomic norm minimization for more general measurement models. While there have been significant advances in the understanding of atomic norm minimization for line spectrum estimation, its theory and application to other measurement settings require further investigation.

• •

Applications in communications, sensing, and imaging. Atomic norm minimization has emerged as a popular approach recently for many practical applications, such as channel estimation in massive MIMO [100, 101], radar imaging [102], and nuclear magnetic resonance (NMR) spectroscopy [103]. It is our hope that this paper will stir more interest in applying the atomic norm in applications that call for high-resolution parameter estimation.

Acknowledgement

The authors thank the editorial board for early feedbacks of the white paper, and the associate editor and anonymous reviewers for constructive suggestions that help improve the quality of this paper. This work is supported in part by ONR under grants N00014-18-1-2142 and N00014-19-1-2404, and by NSF under grants CIF-1826519 and ECCS-1818571.