The dual approach to non-negative super-resolution: impact on primal reconstruction accuracy

TL;DR

This paper investigates the stability of non-negative super-resolution solutions by analyzing the dual approach, establishing a relationship between perturbations in dual and primal variables using advanced mathematical techniques.

Contribution

It provides a novel stability analysis of the dual approach in non-negative super-resolution, linking dual perturbations to primal solution accuracy.

Findings

Established a quantitative relationship between dual and primal perturbations.

Demonstrated stability of super-resolution solutions under dual problem perturbations.

Applied a non-trivial implicit function theorem to analyze solution stability.

Abstract

We study the problem of super-resolution, where we recover the locations and weights of non-negative point sources from a few samples of their convolution with a Gaussian kernel. It has been recently shown that exact recovery is possible by minimising the total variation norm of the measure. An alternative practical approach is to solve its dual. In this paper, we study the stability of solutions with respect to the solutions to the dual problem. In particular, we establish a relationship between perturbations in the dual variable and the primal variables around the optimiser. This is achieved by applying a quantitative version of the implicit function theorem in a non-trivial way.

Equations

$$x=\sum _{k=1}^K{a}_k{\delta }_{t_k},$$

$$y_j={\int }_I\varphi (t-s_j)x(dt),$$

$$y_j=\sum _{k=1}^K{a}_k\varphi (t_k-s_j),$$

$$\underset{\hat{x}\ge 0}{\min }{\int }_I\hat{x}(dt)\text{subject to}\mathbf{y}={\int }_I\Phi (t)\hat{x}(dt).$$

$$\underset{\lambda \in {\mathbb{R}}^M}{\max }{\mathbf{y}}^T\lambda \text{subject to}\sum _{j=1}^M{\lambda }_j\varphi (t-s_j)\le 1,\forall t\in I.$$

$$q(s)=\sum _{j=1}^M{\lambda }_j^{\ast }\varphi (s-s_j),$$

$$q(t_i)=1,\forall i=1,\dots ,K,$$

$$q(s)<1,\forall s\ne t_i,\forall i=1,\dots ,K.$$

$$|t-t^{\ast }|\le C_{t^{\ast }}\Vert \lambda -{\lambda }^{\ast }{\Vert }_2,$$

$${\delta }_{\lambda }\le \frac{|q^{\prime \prime }(t^{\ast }){|}^2{\sigma }^3\sqrt{e}}{4\sqrt{2}\left(2+cR\right)M},$$

$$C_{t^{\ast }}$$

$$R$$

$${\Phi }_{ij}=\varphi (t_j-s_i),$$

$${\mathbf{t}}^{\ast }=[t_1,\dots ,t_K{]}^T,{\mathbf{a}}^{\ast }=[a_1,\dots ,a_K{]}^T.$$

$$\Phi \mathbf{a}=\mathbf{y}.$$

$$\tilde{\Phi }=\Phi +E,$$

$$\underset{\mathbf{a}}{\min }\Vert \tilde{\Phi }\mathbf{a}-\mathbf{y}{\Vert }_2.$$

$$\Vert {\mathbf{t}}^{\ast }-\tilde{\mathbf{t}}{\Vert }_2<\frac{{\sigma }^2{\sigma }_{\max }(\Phi )}{4e^{4/{\sigma }^2}\sqrt{M}}\left(\sqrt{1+\frac{{\sigma }_{\min }^2(\Phi )}{{\sigma }_{\max }^2(\Phi )}}-1\right).$$

$$\Vert {\mathbf{a}}^{\ast }-\tilde{\mathbf{a}}{\Vert }_2\le \frac{4e^{\frac{4}{{\sigma }^2}}\sqrt{M}\Vert {\mathbf{a}}^{\ast }{\Vert }_2}{{\sigma }^2{\sigma }_{\min }(\Phi )}\Vert \tilde{\mathbf{t}}-{\mathbf{t}}^{\ast }{\Vert }_2+O(\Vert \tilde{\mathbf{t}}-{\mathbf{t}}^{\ast }{\Vert }_2^2).$$

$$F(t,\lambda )=\sum _{j=1}^M{\lambda }_j{\varphi }^{\prime }(t-s_j),$$

$${\partial }_{\lambda }t(\lambda )=-{\left[{\partial }_{t}F(t(\lambda ),\lambda )\right]}^{-1}{\partial }_{\lambda }F(t(\lambda ),\lambda ),$$

$$\underset{(t,\lambda )\in V_{\delta }}{\sup }\Vert I-{\left[{\partial }_{t}F(t^{\ast },{\lambda }^{\ast })\right]}^{-1}{\partial }_{t}F(t,\lambda )\Vert \le \frac{1}{2},$$

$${\delta }_1=(2M_t{B}_{\lambda }{)}^{-1}{\delta }_0,$$

$$B_{\lambda }$$

$$M_t$$

$${\delta }_0=\frac{{\sigma }^2|q^{\prime \prime }(t^{\ast })|}{\sqrt{M}\left(4+2c\cdot \frac{\Vert {\lambda }^{\ast }{\Vert }_2}{\sigma }\right)}.$$

$${\delta }_{\lambda }=\frac{\sigma \sqrt{e}|q^{\prime \prime }(t^{\ast })|}{2\sqrt{2M}}\cdot {\delta }_0.$$

$${\partial }_{t}F(t,\lambda )$$

$${\partial }_{\lambda }F(t,\lambda )$$

$$t(\lambda )=t({\lambda }^{\ast })+⟨\lambda -{\lambda }^{\ast },{\partial }_{\lambda }t({\lambda }_{\delta })⟩,$$

$$\left|t(\lambda )-t({\lambda }^{\ast })\right|\le {\Vert \lambda -{\lambda }^{\ast }\Vert }_2\cdot {\Vert {\partial }_{\lambda }t({\lambda }_{\delta })\Vert }_2.$$

$${\delta }_t\le \frac{C+{\delta }_{t}A}{B}\cdot {\delta }_{\lambda }.$$

Full text

I Problem setup

In the study of non-negative super-resolution, we want to estimate a signal $x$ which consists of a number of point sources with unknown locations and non-negative magnitudes, from only a few measurements of the convolution of $x$ with a known kernel $\phi$. This is a problem that arises in a number of applications, for example fluorescence microscopy [8], astronomy [9] and ultrasound imaging [10]. In such applications, the measurement device has limited resolution and cannot distinguish between distinct point sources in the input signal $x$. This is often modelled as a deconvolution problem with a Gaussian kernel.

Specifically, let $x$ be a non-negative measure on $I=[0,1]$ consisting of $K$ unknown non-negative point sources:

[TABLE]

with $a_{k}>0$, for all $k=1,\ldots,K$, and let $\mathbf{y}\in\mathbb{R}^{M}$ be the vector of measurements obtained by sampling the convolution of $x$ with a known kernel $\phi$ (e.g. a Gaussian) at locations $s_{j}$:

[TABLE]

for all $j=1,\ldots,M$. Note that, because $x$ is a discrete measure, each entry in $\mathbf{y}$ is of the form:

[TABLE]

for all $j=1,\ldots,M$. Let $\Phi(t)=[\phi(t-s_{1}),\ldots,\phi(t-s_{M})]^{T},$ then $x$ can be recovered by solving the following program: 111 We assume that the measurements are exact. We treat the case when the measurements are noisy in the journal version of this paper.

[TABLE]

The problem of super-resolution has been studied extensively in the literature since the seminal paper [1], which addressed the case of complex amplitudes. In [2], the authors showed that when the coefficients $a_{k}$, $k=1,\ldots,K$ are positive, which corresponds to the present setting, exact recovery is possible without separation of sources.

The dual of problem (1) is

[TABLE]

The dual problem (2) is a finite-dimensional problem with infinitely many constraints, known as a semi-infinite program. Such problems can be solved using a number of algorithms including exchange methods [7] and sequential quadratic programming [6]. The advantage over algorithms that solve the primal problem (for example the ADCG algorithm [3]) is working in a finite dimensional space, which simplifies the analysis.

Consider a solution $\lambda^{*}$ of the dual problem (2) which corresponds to a dual certificate, namely a function

[TABLE]

which satisfies the conditions:

[TABLE]

Then the local maximisers of $q(s)$ correspond to the source locations $\{t_{k}\}_{k=1}^{K}$ and the amplitudes $\{a_{k}\}_{k=1}^{K}$ are found by solving a linear system.

In this paper, we analyse how small perturbations of $\lambda^{*}$ affect the local maximisers of $q(s)$ in the case when the convolution kernel is Gaussian $\phi(t)=e^{-t^{2}/\sigma^{2}}$. The outcome is a bound on how far the estimated locations $t_{k}$ and magnitudes $a_{k}$ are from their true values obtained for exact $\lambda^{*}$. This gives us an insight into the size of the error in the locations and magnitudes when we apply an optimisation algorithm to the dual of the super-resolution problem.

II Bound on the error as $\lambda$ is perturbed

In this section we present our main results, namely two theorems that give bounds on the perturbations around the source locations $t_{k}$ and the magnitudes $a_{k}$ respectively, as the dual variable is perturbed away from the optimiser $\lambda^{*}$, when the convolution kernel is a Gaussian with known width $\sigma$: $\phi(t)=e^{-t^{2}/\sigma^{2}}.$

Theorem 1

(Dependence of $t$ on $\lambda$) Let $\lambda^{*}\in\mathbb{R}^{M}$ be a solution of the dual program (2) with $\phi$ Gaussian as given above such that the dual certificate $q(s)$ defined in (3) satisfies conditions (4) and (5), $\lambda$ a perturbation of $\lambda^{*}$ in a ball of radius $\delta_{\lambda}$ and $t$ an arbitrary local maximiser of $q_{\lambda}(s)=\sum_{j=1}^{M}\lambda_{j}\phi(s-s_{j})$ so that for $\lambda=\lambda^{*}$, the corresponding local maximiser $t^{*}$ is a source location of $x$. Then

[TABLE]

provided that

[TABLE]

where

[TABLE]

and $c\approx 3.9036$ is a universal constant.

One of the main conclusions which can be drawn from this result is that the primal spike location error is controlled in $l_{\infty}$, but degrades as a function of the number of measurements in the order of $\sqrt{M}$. Of crucial importance is the curvature of the dual certificate at the true solution: the flatter the certificate, the worse the estimation error. Our theorem also gives important information about the accuracy in the dual variable required to guarantee our upper bound on the error of recovery. This accuracy is of the inverse order of the number of measurements, which is quite a stringent constraint. Both the $M$ and the $\sqrt{M}$ factors are a consequence of the way we bound sums of shifted copies of the kernel, namely $\sum_{j=1}^{M}\phi(t-s_{j})\leq M\max_{t\in\mathbb{R}}\phi(t)$. Given the fast decay of the Gaussian, it is clear that this is not a tight bound. However, any bound would reflect the density of samples close to each source location.

We will now give a result regarding the perturbation of the magnitudes $a_{k}$ when $\lambda^{*}$ is perturbed. Let $\Phi$ be the matrix whose entries are defined as

[TABLE]

and $\mathbf{t^{*}}$ and $\mathbf{a^{*}}$ the vectors of true source locations and weights:

[TABLE]

When we solve (2) exactly, we obtain the source locations by finding the local maximisers of $q(s)$. Then the vector of weights $\mathbf{a^{*}}$ is found by solving the system

[TABLE]

When the source locations are perturbed, we denote the resulting perturbed data matrix by:

[TABLE]

and we calculate the vector of perturbed weights $\tilde{\mathbf{a}}$ as the solution of the least squares problem

[TABLE]

The following theorem gives us a bound on the difference $\|\mathbf{a^{*}}-\tilde{\mathbf{a}}\|_{2}$ between the vector of true weights $\mathbf{a^{*}}$ and the vector of weights $\tilde{\mathbf{a}}$ obtained by solving the least squares problem (10) with the perturbed matrix $\tilde{\Phi}$, as a function of the difference $\|\tilde{\mathbf{t}}-\mathbf{t^{*}}\|_{2}$ between the perturbed source locations $\tilde{\mathbf{t}}$ and the true source locations $\mathbf{t^{*}}$.

Theorem 2

(Dependence of $a$ on $t$) Let $\mathbf{t^{*}}\in[0,1]^{K}$ be the vector of true source locations and $\tilde{\mathbf{t}}\in[0,1]^{K}$ the perturbed source locations, such that:

[TABLE]

Then the error between the true weights $\mathbf{a^{*}}$ and the perturbed weights $\tilde{\mathbf{a}}$ obtained by solving program (10) is bounded by:

[TABLE]

III Proofs

In this section we present the proofs of the two theorems. Due to space limitations, we skip some details, which will be present in the journal version of the paper.

III-A Proof of Theorem 1

Let $t^{*}$ be an arbitrary local maximiser of the function $q(t)$ in (3), so $t^{*}$ is also a source location, and $\lambda^{*}$ the solution to (2). The key step in this proof is applying a quantitative version of the Implicit Function Theorem [11] to the function:

[TABLE]

where $F(t^{*},\lambda^{*})=0$ because $t^{*}$ is a maximizer of $q(s)$ in (3). The theorem tells us that we can express $t$ as a function $t(\lambda)$ of $\lambda$ with:

[TABLE]

for $t$ in a ball of radius $\delta_{0}$ around $t^{*}$ and for $\lambda$ in a ball of radius $\delta_{1}\leq\delta_{0}$ around $\lambda^{*}$, where $\delta_{0}$ is chosen such that

[TABLE]

where $V_{\delta}=\left\{(t,\lambda)\in\mathbb{R}^{M+1}\mathrel{\mathop{\mathchar 58\relax}}|t-t^{*}|\leq\delta_{0},\|\lambda-\lambda^{*}\|\leq\delta_{0}\right\}$ and $\delta_{1}$ is given by

[TABLE]

where

[TABLE]

The following two lemmas give values of $\delta_{0}$ and $\delta_{1}$ that define balls around $t^{*}$ and $\lambda^{*}$ respectively which are included in the balls required by the Quantitative Implicit Function Theorem with radii defined in (14) and (15).

Lemma 1

(Radius of ball around $t^{*}$) The condition (14) is satisfied if

[TABLE]

Lemma 2

(Radius of ball around $\lambda^{*}$) For $\delta_{0}$ from Lemma 1 and $\delta_{1}$ from condition (15), we have that $\delta_{\lambda}<\delta_{1}$ if

[TABLE]

Given the definition of the function $F$ in (12), we have that

[TABLE]

By applying a Taylor expansion to $t(\lambda)$ around $\lambda^{*}$ in the region defined by $\delta_{0}$ and $\delta_{\lambda}$, we have that

[TABLE]

for some $\lambda_{\delta}$ on the line segment determined by $\lambda^{*}$ and $\lambda$, so

[TABLE]

Bounding the second factor in (16) above, using (13) and applying a number of manipulations, we obtain

[TABLE]

where $\delta_{t}=|t(\lambda)-t(\lambda^{*})|$, $\|\lambda-\lambda^{*}\|_{2}\leq\delta_{\lambda}$ and

[TABLE]

We now need to lower bound $B$ and upper bound $C+\delta_{t}A$. For $B$, by applying a Taylor expansion around $t(\lambda^{*})-s_{j}$, each term in the sum is equal to

[TABLE]

for $j=1,\ldots,M$ and some $\xi_{j}$ in the interval

[TABLE]

By combining this with the reverse triangle inequality, Cauchy-Schwartz inequality and Lemma 1, we obtain

[TABLE]

where $c\approx 3.9036$. By using the global maximum of $\phi^{\prime}(t)$ and $\phi^{\prime\prime}(t)$ for $\phi(t)$ Gaussian, we have that

[TABLE]

By finally plugging the above inequalities and the result of Lemma 1 into (17), we obtain the expression for $C_{t^{*}}$ in (6).

III-B Proof of Theorem 2

We apply equation (4.2) from [5] with $e=0$ (the noise in the observations) and obtain

[TABLE]

where $\Phi^{\dagger}=(\Phi^{T}\Phi)^{-1}\Phi^{T}$ is the pseudoinverse of $\Phi$ and $F=O(E)$ is the perturbation of the $\Phi^{\dagger}$ due to the perturbation $E$ of $\Phi$, namely

[TABLE]

In order to obtain an explicit expression for $F$, we write $\tilde{\Phi}^{\dagger}$:

[TABLE]

where

[TABLE]

We then write the first factor in (19) as

[TABLE]

where

[TABLE]

and in the second inequality in (21) we applied the Neumann series expansion to the matrix $-{\Phi^{\dagger}}^{T}\Delta\Phi^{\dagger}$, which converges if

[TABLE]

We will return to condition (23) at the end of this section. We now substitute (21) in (19), giving

[TABLE]

so we have that

[TABLE]

which is indeed $O(E)$, since $S_{\Phi}=O(\Delta)$ and $\Delta=O(E)$. We next upper bound $\|S_{\Phi}\|_{2}$. From (22) we have

[TABLE]

Now let $D$ be an upper bound on $\|\Delta\|_{2}$, obtained by applying the triangle inequality in (20), so that

[TABLE]

Then, from (25) we have

[TABLE]

where the series converges if $D\|\Phi^{\dagger}\|_{2}^{2}<1$, in which case the denominator in the last fraction above is positive. We return to this condition at the end of the section. We also know that222 Using the SVD $\Phi=U\Sigma V^{T}$, we have $\Phi^{\dagger}=(\Phi^{T}\Phi)^{-1}\Phi^{T}=(V\Sigma^{2}V^{T})^{-1}V\Sigma U^{T}=V\Sigma^{-1}U^{T}$, so the conclusion follows.

[TABLE]

By applying the triangle inequality in (24) and then using (27) and the fact that $\|(\Phi^{T}\Phi)^{-1}\|_{2}=1/\sigma_{\min}^{2}(\Phi)=\|\Phi^{\dagger}\|_{2}^{2}$ (from (28)), we obtain

[TABLE]

where $D$ is given in (26). It remains to establish an upper bound on $\|E\|_{F}$, and consequently on $\|E\|_{2}$. The following lemma gives us such a bound.

Lemma 3

(Upper bound of $\|E\|_{F})$ Let $E=\tilde{\Phi}-\Phi$ for $\Phi$ and $\tilde{\Phi}$ as defined in (8) and (9) respectively for $t_{j},\tilde{t}_{j}\in[0,1]$ for $j=1,\ldots,K$. Then:

[TABLE]

By using the triangle inequality and norm sub-multiplicativity in (18), and then substituting (29) and (30), we obtain

[TABLE]

which is the bound given in Theorem 2. Note that because $\|E\|_{2}=O(\|\tilde{\mathbf{t}}-\mathbf{t^{*}}\|_{2})$ (see (30)), the first term is the only term that is $O(\|\tilde{\mathbf{t}}-\mathbf{t^{*}}\|_{2})$ in the first inequality above, so the other terms are included in the $O(\|\tilde{\mathbf{t}}-\mathbf{t^{*}}\|_{2}^{2})$ term at the end.

Lastly, we return to condition (23), which must be satisfied in order for the bound above to hold. By using norm sub-multiplicativity and the bound on $\|\Delta\|_{2}$ from (26), we obtain

[TABLE]

and by requiring that the right hand side above is less than one, we obtain a quadratic constraint on $\|E\|_{2}$, satisfied if

[TABLE]

By using the bound on $\|E\|_{2}$ from (30), the above holds if (11) holds. Note that by imposing this, we also ensure that the condition for the series in (27) to converge holds, since $D\|\Phi^{\dagger}\|_{2}^{2}$ is equal to the right hand side of (31).

IV Conclusion

In this paper, we proved primal stability in the non-negative super-resolution problem, when addressed via convex duality. The main ingredient in our analysis is a quantitative version of the implicit function theorem, a folklore result in the theory of dynamical systems community.

Our results provide precise orders in the number of measurements for the accuracy of the solution to the convex dual problem and an $\ell_{\infty}$ error bound on the primal spike locations.

Future plans include the study of the dual approach to the noisy super-resolution problem using similar techniques to the ones developed here.