arXiv:1501.01825v1 [cs.IT] 8 Jan 2015 Unified Convex Optimization Approach to Super-Resolution Based on Localized Kernels Tamir Bendory Shai Dekel Electrical Engineering Technion - Israel Institute of Technology GE Global Research School of Mathematical Sciences Tel-Aviv University Abstract—The problem of resolving the fine details of a signal from its coarse scale measurements or, as it is commonly referred to in the literature, the super-resolution problem arises naturally in engineering and physics in a variety of settings. We suggest a unified convex optimization approach for super-resolution. The key is the construction of an interpolating polynomial based on localized kernels. We also show that the localized kernels act as the connecting thread to another wide-spread problem of stream of pulses. I. I NTRODUCTION Super-resolution is the task of estimating the fine details of a signal from its coarse scale measurements (see for instance [1], [2], [3]). This problem arises in many physical and engineering problems as sensing systems have a physical constraint on the highest resolution the system can achieve. The aim of this paper is to suggest a unified convex optimization approach for super-resolution problems, based on the construction of interpolating polynomials and the existence of well-localized kernels. For simplicity, we consider signals of the form X x(t) = cm δtm , t ∈ M, (I.1) m where δt is a Dirac measure, M ⊂ Rd , d ≥ 1, is a compact manifold and T := {tm } is the signal’s support. This model can be extended to piece-wise polynomials given mild additional data (e.g. boundary condition, signal’s average). The information we have on the signal is its low spectral coefficients, namely its ’low-resolution’ measurements. In this paper, we consider three concrete examples. Suppose M is the d-dimensional torus. In this case, we assume that the measurements are given as y = FN x, where FN projects x onto the space of trigonometric polynomials of degree N . That is, we have access solely the low 2N + 1 low Fourier coefficients of x. This model reflects the fact that many measurement devices have limited resolution (e.g. diffraction limits in microscopes) . This problem can be solved using parametric methods such as MUSIC, matrix pencil and ESPRIT [4], [5], [6], [7]. However, these methods tend to be unstable in the presence of noise or model mismatch due to sensitivity of polynomial root finding. Arie Feuer Electrical Engineering Technion - Israel Institute of Technology Another closely-related example is the projection onto the space of algebraic polynomials. In this case we choose M as the manifold [−1, 1]d . This problem arises in spectral methods to numerically solve partial differential equations (see e.g. [8], [9] ). The third example is low-resolution measurements of signals lying on the bivariate sphere S2 . In this case, we assume that a Dirac ensemble of the form (I.1) lies on S2 , and we have access only to the low N spherical harmonics, which is the extension of Fourier analysis for signals on the sphere. This problem arises for instance in medical imaging [10], [11]. The rest of the paper is organized as follows. The following section describes the convex optimization framework and its main ingredient - the construction of interpolating polynomials. Section III reveals the intimate relations between localized kernels and super-resolution and demonstrates it on the three examples. Later on, Section IV is devoted to the problem of recovery from stream of pulses and its connection to the super-resolution problem. Section V shows some numerical experiments for super-resolution on the sphere. Ultimately, in Section VI we draw some conclusions and suggest future extensions. II. C ONVEX O PTIMIZATION A PPROACH TO S UPER - RESOLUTION In this paper we focus on a convex optimization approach for resolving signals from their low resolution measurements. We use the Total-Variation (TV) norm as a sparse-promoting regularization. In essence, the TV norm is the generalization of `1 norm to the real line (for rigorous definition, see for instance Section 1 in [12]Pand [13]). For signals of the form (I.1), we have kxkT V = m |cm |. The main pillar of this framework is the following theorem [14]: P Theorem II.1. Let x(t) = m cm δtm , cm ∈ R, where T := {tm } ⊆ M , and M is a compact manifold in Rn . Let ΠD be a linear space of continuous functions of dimension D in M . For any basis {Pk }D k=1 of ΠD , let yk = hf, Pk i for all 1 ≤ k ≤ D. If for any set {um }, um ∈ R, with |um | = 1, there exists q ∈ ΠD such that q(tm ) = um , ∀tm ∈ T, (II.1) |q(t)| < 1 , ∀t ∈ M \T, (II.2) then x is the unique real Borel measure satisfying min x ˜∈M(M ) k˜ xkT V subject to yk = h˜ x, Pk i , 1 ≤ k ≤ D. (II.3) By taking ΠD as the space spanned by the eigenfunctions of the low-resolution measurement operator, Theorem II.1 implies that the super-resolution problem can be reduced to a polynomial construction problem. We mention that similar techniques, frequently referred to as dual certificate, are widely used in many sparse recovery problems, see for instance [15], [16], [17]. From the algorithmic perspective, it is interesting to notice that in all three cases the infinite-dimensional TV minimization problem (II.3) can be solved accurately (i.e. to any desired resolution) by algorithms with finite complexity which depends only on N . The recovery is performed by a three-stage algorithm involving the solution of a semi-definite program based on the dual problem of (II.3) (can be solved using offthe-shelf software), followed by root-finding and least square fitting [18], [12], [19]. Alternatively, one can discretize the problem and solve standard `1 minimization which converges to the solution of (II.3) as the discretization becomes finer [20]. As we discuss later, these solutions are robust to noisy measurements under the separation condition. III. T HE C ONSTRUCTION OF I NTERPOLATING P OLYNOMIALS As aforementioned, super-resolution problem can be reduced to the construction of an interpolating polynomial which lies in the space spanned by the eigenfunctions of the lowresolution measurement operator. The existence of such a polynomial relies on two interrelated pillars. The first is the separation condition defined as follows. Definition III.1. A set of points T ⊂ M is said to satisfy the minimal separation condition with respect to the metric d(·, ·) if ∆ := min d (ti , tj ) ≥ ν/N, ti ,tj ∈T,ti 6=tj where ν > 0 is a constant which does not depend on N . Along this work we argue that the separation condition is a sufficient condition. However, we emphasize that few works prove that the separation is also necessary, and without minimal separation the recovery can not be stable in the presence of superfluous noise [18], [19], [21]. The second pillar is the existence of well-localized kernels in the space spanned by the eigenfunction of the low-resolution measurement operator. To explain this statement, we will demonstrate it through our three prime examples. Consider a signal of the form of (I.1) defined on the circle [0, 1] and its projection onto the space of trigonometric polynomials of degree N . In the time domain, the projection is given as y(t) = (x ∗ DN )(t), PN where DN (t) = k=−N ejkt is the Dirichlet kernel of degree N. The key for super-resolving the signal is the existence ˜ N , a well-localized smooth super-position the Fejer kernel D of Dirichlet kernels. Equipped with the wrap-around metric d (ti , tj ) = kti − tj k∞ , the authors of [18] showed that under the separation condition of Definition III.1, there exist coefficients {am } and {bm } so that the polynomial X ˜ N (t − tm ) + bm D ˜ (1) (t − tm ) , q(t) = am D N m satisfies (II.1) and (II.2). Hence, the TV minimization (II.3) recovers x exactly in the univariate case. They also showed that similar polynomials can be reconstructed for higher dimensional signals. In this manner, the existence of the Fejer kernel is crucial for the ability to resolve a signal on the d-dimensional torus under the separation condition. In consecutive papers, the existence of the interpolating polynomial also serve as the basis for proving the robustness and localization of the solution of the TV minimization (II.3) to noisy measurements [22], [23]. A similar result holds for the case in which the measurements are the projection of a signal of the form (I.1) on [−1, 1] onto any basis spanning the space of algebraic polynomials of degree N . As was shown in [12], by considering the metric d (ti , tj ) = |arccos (ti ) − arccos (tj )| one can construct the appropriate algebraic polynomial q(t) under the separation condition of Definition III.1 and hence the recovery through TV minimization is guaranteed. Observe that in this case, the minimal required separation is space-dependent and reduced near the edges to the order of O N −3/2 . A consecutive paper showed that the recovery is also robust to noisy measurements [24]. These results hold for the bivariate case as well. The same phenomenon occurs for signals on the bivariate sphere S2 . Recall that spherical harmonics is the natural extension of Fourier analysis to signals on the sphere, and consider a Dirac ensemble on S2 X x(ξ) = cm δξm , {ξm } ∈ Ξ ⊂ S2 . (III.1) m Let Yn,k (ξ), 0 ≤ n ≤ N, −n ≤ k ≤ n be an orthonormal basis of VN S2 , the space of spherical harmonics of degree N . We assume that the information we have on the signal is its projection onto VN S2 yˆn,k = hx, Yn,k i, 0 ≤ n ≤ N, −n ≤ k ≤ n. (III.2) In the space domain, the projection onto VN S2 can be written as a spherical convolution, Z y(ξ) = x(η)KN (ξ · η)dS2 (η), S2 where by the addition formula [25] KN (ξ · η) = N X k=−N Yn,k (ξ)Y n,k (η) = 2N + 1 PN,3 (ξ · η), 4π where Y n,k is the conjugate of Yn,k , and PN,3 (x) is the univariate ultraspherical Gegenbauer polynomial of order 3 and degree N . In [26], it was shown that a smooth super-position of the kernel KN (ξ · η) results in a well-localized kernel in the space of spherical harmonics of degree N . We denote the localized kernel as FN (ξ ·η), and by Dξ,` , ` = 1, 2 the partial rotational derivatives at ξ. Leveraging the localization of FN , it is known [14] that there exist coefficients {αm }, {βm } and {γm } so that a polynomial of the form X q(ξ) = αm FN (ξ · ξm ) + βm Dξm ,1 FN (ξ, ξm ) m + γm Dξm ,2 FN (ξ, ξm ), fulfils the requirements of Theorem II.1 under the separation condition with the natural distance on the sphere d (ξi , ξj ) = arccos (ξi · ξj ). Accordingly, we present the following theorem: Theorem III.2. [14] Let Ξ = {ξm } be the support of a signed measure x of the form (III.1). Let {Yn,k }N n=0 be any spherical harmonics basis for VN (S2 ) and let yˆn,k = hx, Yn,k i, 0 ≤ n ≤ N , −n ≤ k ≤ n. If Ξ satisfies the separation condition of Definition III.1 with the metric d (ξi , ξj ) = arccos (ξi · ξj ), then x is the unique solution of min x ˜∈M(S2 ) k˜ xkT V problem of stream of pulses, namely recovery of the delays T := {tm } and weights {cm } from the measurements X y(t) = cm Kσ (t − tm ) , cm , t ∈ R, m where K(t) is a pulse shape (kernel) and Kσ (t) := K(t/σ) for some σ > 0. An alternative representation of the problem is as y(t) = (x ∗ Kσ ) (t) where X x(t) = cm δtm , T := {tm } , m and δt is a Dirac measure. In a manner similar to Theorem II.1, the recovery problem can be reduced to construction of a special interpolating function, as follows: P Theorem IV.1. [27] Let x(t) = R m cm δtm , cm ∈ R where T := {tm } ⊆ R, and let y(t) = R K(t − s)dx (s) for a L times differentiable kernel K(t). If for any set {um }, um ∈ R, with |um | = 1, there exists a function of the form Z X L q(t) = K (`) (s − t)dµ` (s) , R `=0 L for some measures {µ` (t)}`=0 , satisfying q(tm ) = um , ∀tm ∈ T, |q(t)| < 1 , ∀t ∈ R\T, subject to yˆn,k = h˜ x, Yn,k i (III.3) 0 ≤ n ≤ N, −n ≤ k ≤ n, then x is the unique real Borel measure solving min k˜ xkT V subject to Z y(t) = K(t − s)d˜ x (s) . where M(S2 ) is the space of signed Borel measures on S2 . It has been shown in [19] (for a discrete version of (III.3)) that the recovery error in the presence of noise is proportional to the noise level. We conclude this section with an interesting observation regrading non-negative signals (i.e. cm > 0). In this case, a sufficient condition for signal recovery is the existence of interpolating polynomial as in Theorem II.1 but the constraint in (II.1) is replaced by a weaker constraint of q(tm ) = 1 for all tm ∈ T (see Theorem 5.1 in [14]). Consider the case of the projection signals on [0, 1] onto the space of trigonometric polynomials of degree N . In this case for any s ≤ N a polynomial of the form q(t) = 1 − 2−(s+1) s Y (1 − cos (t − tm )) , m=1 is a trigonometric polynomial of degree N , q(tm ) = 1 for all tm ∈ T and |q(t)| < 1 otherwise. Consequently, a clustered N -sparse signal with non-negative coefficient can be recovered exactly. The same observation is noted for the other two cases with similar constructions (see [14]). IV. S UPER - RESOLUTION AND S TREAM OF P ULSES In the previous section we stressed that the existence of well-localized kernels is crucial for resolving signals. The localized kernels bind the super-resolution problem to the x ˜∈M(R) (IV.1) R We note that here too, as in the super-resolution problems, the existence of q(t) relies on separation and localization. If the support T satisfies the separation condition with the metric d(ti , tj ) = |ti − tj | and σ = 1/N , there exists coefficients {am } and {bm } so that X q(t) = am Kσ (t − tm ) + bm Kσ(1) (t − tm ) , m satisfies the requirements of Theorem IV.1, and thus enabling perfect recovery. This holds if the kernel K(t) is localized. More precisely, the kernel should satisfy the definition of admissible kernel as follows: Definition IV.2. A kernel K is admissible if it has the following properties: 1) K ∈ C 3 (R), is real and even. 2) Global property: There (`) exist constants C`2 > 0, ` = 0, 1, 2, 3 such that K (t) ≤ C` / 1 + t . 3) Local property: There exist constants ε, β > 0 such that a) K(t) > 0 for all |t| ≤ ε and K(t) < K(ε) for all |t| > ε, b) K (2) (t) < −β for all |t| ≤ ε. Here too, the existence of the interpolating function guarantees a robust and localized recovery [27], [28]. Fig. 1: The mean recoery error (in log scale) as a function of ν over 20 simulations. To be clear, by error we merely mean the distance on the sphere between the true and the estimated support. (a) The low resolution measurements. V. N UMERICAL E XPERIMENTS ON THE S PHERE To demonstrate our results we chose to use the less familiar example where we consider the recovery of a Dirac ensemble on the sphere (III.1) from its low-resolution measurements (III.2). In this case, the TV minimization (III.3) can be solved to any desired accuracy by a three-stages algorithm consists of semi-definite optimization problem with O N 4 variables, followed by root-finding and least square fitting. In case that the measurements are contaminated with bounded noise, it has been shown that the recovery error is proportional to the noise level (in the discrete setting) [19]. The following experiments were conducted in Matlab using CVX [29], which is the standard modelling system for convex optimization. The signals were generated in the following two stages: • • Random locations on the sphere were drawn uniformly, sequentially added to the signal’s support, while maintaining the separation condition of Definition III.1. Once the support was determined, the amplitudes were drawn randomly from an iid normal distribution with standard deviation of SD = 10. Figure 1 presents a numerical estimation of the separation constant ν (see Definition III.1). As can be seen, a separation constant of 2π seems to be sufficient in the noise-free setting. This separation coincides with the spatial resolution of the projection of x onto VN [30]. Figure 2 presents an example for super-resolution on the sphere. The signal in Figure 2a is the projection of the signal onto V10 , whereas Figure 2b presents the recovered signal. The maximal and average recovery errors for several values of N are presented in Table I. In the noisy setting, we added an iid noise with normal distribution and zero mean. Figure 3 presents the mean recovery error as a function of the noise standard deviation. We note that the error increases moderately with the standard deviation. (b) The recovered signal f . Fig. 2: Super-resolution on the sphere for N = 10. The signal is presented on a grid for visualization only. N Average error Max error 5 8.1267 × 10−5 2.163 × 10−4 8 8.1826 × 10−5 1.9 × 10−3 10 9.0404 × 10−5 3.3 × 10−3 TABLE I: The localization error for N = 5, 8, 10. For each value of N , 10 experiments were conducted. To be clear, by error we merely mean the distance on the sphere between the true and the estimated support. VI. D ISCUSSION In this paper, we have presented a general framework for resolving robustly signals in various settings and geometries using convex optimization. Localized kernels play a crucial role in the process. An important question is how general is this framework. That is to say, under which settings we can expect to resolve signal robustly using convex optimization. This topic is under ongoing research. We also showed that the localization principle relates superresolution to other problems such as the recovery from stream of pulses. Recently, some new results on the localization of the solution of (IV.1) have been developed [28] and we plan to test it on real ultrasound data. It is of a great interest to find Fig. 3: 10 experiments were conducted with N = 8 for each value of standard deviation. The figure presents the average recovery error as a function of the noise standard deviation. By error, we merely mean the distance on the sphere between the true and the estimated supports. more applications where similar techniques could be applied (for recent related works, see for instance [31], [32]). From the algorithmic perspective, the aforementioned superresolution problems can be recast as finite dimensional problems involve the solution of a semi-definite program, rootfinding and least square fitting. It will be interesting to look for the relations of these algorithms with ’traditional’ parametric methods such as MUSIC which also rely on root-finding. R EFERENCES [1] H. Greenspan, “Super-resolution in medical imaging,” The Computer Journal, vol. 52, no. 1, pp. 43–63, 2009. [2] S. C. Park, M. K. Park, and M. G. Kang, “Super-resolution image reconstruction: a technical overview,” Signal Processing Magazine, IEEE, vol. 20, no. 3, pp. 21–36, 2003. [3] M. Elad and A. Feuer, “Restoration of a single superresolution image from several blurred, noisy, and undersampled measured images,” Image Processing, IEEE Transactions on, vol. 6, no. 12, pp. 1646–1658, 1997. [4] P. Stoica and R. L. Moses, Spectral analysis of signals. Pearson/Prentice Hall Upper Saddle River, NJ, 2005. [5] Y. Hua and T. Sarkar, “Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise,” Acoustics, Speech and Signal Processing, IEEE Transactions on, vol. 38, no. 5, pp. 814–824, 1990. [6] R. Roy and T. Kailath, “Esprit-estimation of signal parameters via rotational invariance techniques,” Acoustics, Speech and Signal Processing, IEEE Transactions on, vol. 37, no. 7, pp. 984–995, 1989. [7] R. Schmidt, “Multiple emitter location and signal parameter estimation,” Antennas and Propagation, IEEE Transactions on, vol. 34, no. 3, pp. 276–280, 1986. [8] C. Canuto, M. Hussaini, A. Quarteroni, and T. Zang, Spectral Methods in Fluid Dynamics (Scientific Computation). Springer-Verlag, New York-Heidelberg-Berlin, 1987. [9] J. Shen, “Efficient spectral-galerkin method i. direct solvers of secondand fourth-order equations using legendre polynomials,” SIAM Journal on Scientific Computing, vol. 15, no. 6, pp. 1489–1505, 1994. [10] J. Tournier, F. Calamante, D. G. Gadian, A. Connelly et al., “Direct estimation of the fiber orientation density function from diffusionweighted mri data using spherical deconvolution,” NeuroImage, vol. 23, no. 3, pp. 1176–1185, 2004. [11] S. Deslauriers-Gauthier and P. Marziliano, “Spherical finite rate of innovation theory for the recovery of fiber orientations,” in Engineering in Medicine and Biology Society (EMBC), 2012 Annual International Conference of the IEEE. IEEE, 2012, pp. 2294–2297. [12] T. Bendory, S. Dekel, and A. Feuer, “Exact recovery of non-uniform splines from the projection onto spaces of algebraic polynomials,” Journal of Approximation Theory, vol. 182, no. 0, pp. 7 – 17, 2014. [13] W. Rudin, Real and complex analysis (3rd). New York: McGraw-Hill Inc, 1986. [14] T. Bendory, S. Dekel, and A. Feuer, “Exact recovery of dirac ensembles from the projection onto spaces of spherical harmonics,” Constructive Approximation, to appear, 2014. [15] E. Candes, “Mathematics of sparsity (and a few other things),” ICM 2014 Proceedings, to appear, 2014. [16] G. Tang, B. N. Bhaskar, P. Shah, and B. Recht, “Compressive sensing off the grid,” in Communication, Control, and Computing (Allerton), 2012 50th Annual Allerton Conference on. IEEE, 2012, pp. 778–785. [17] Y. De Castro and F. Gamboa, “Exact reconstruction using beurling minimal extrapolation,” Journal of Mathematical Analysis and applications, vol. 395, no. 1, pp. 336–354, 2012. [18] E. J. Cand`es and C. Fernandez-Granda, “Towards a mathematical theory of super-resolution,” Communications on Pure and Applied Mathematics, 2013. [19] T. Bendory, S. Dekel, and A. Feuer, “Super-resolution on the sphere using convex optimization,” submitted. [20] G. Tang, B. N. Bhaskar, and B. Recht, “Sparse recovery over continuous dictionaries-just discretize,” in Signals, Systems and Computers, 2013 Asilomar Conference on. IEEE, 2013, pp. 1043–1047. [21] A. Moitra, “The threshold for super-resolution via extremal functions,” arXiv preprint arXiv:1408.1681, 2014. [22] E. J. Candes and C. Fernandez-Granda, “Super-resolution from noisy data,” Journal of Fourier Analysis and Applications, vol. 19, no. 6, pp. 1229–1254, 2013. [23] C. Fernandez-Granda, “Support detection in super-resolution,” arXiv preprint arXiv:1302.3921, 2013. [24] Y. De Castro and G. Mijoule, “Non-uniform spline recovery from small degree polynomial approximation,” arXiv preprint arXiv:1402.5662, 2014. [25] K. Atkinson and W. Han, Spherical harmonics and approximations on the unit sphere: an introduction. Springer, 2012, vol. 2044. [26] F. Narcowich, P. Petrushev, and J. Ward, “Decomposition of besov and triebel–lizorkin spaces on the sphere,” Journal of Functional Analysis, vol. 238, no. 2, pp. 530–564, 2006. [27] T. Bendory, S. Dekel, and A. Feuer, “Robust recovery of stream of pulses using convex optimization,” submitted. [28] T. Bendory, A. Bar-Zion, S. Dekel, A. Feuer, and D. Adam, “Robust and localized recovery of stream of pulses with application to ultrasound imaging,” in preparation. [29] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” Mar. 2014. [30] B. Rafaely, “Plane-wave decomposition of the sound field on a sphere by spherical convolution,” The Journal of the Acoustical Society of America, vol. 116, no. 4, pp. 2149–2157, 2004. [31] R. Heckel, V. I. Morgenshtern, and M. Soltanolkotabi, “Super-resolution radar,” arXiv preprint arXiv:1411.6272, 2014. [32] C. Aubel, D. Stotz, and H. B¨olcskei, “Super-resolution from short-time fourier transform measurements,” arXiv preprint arXiv:1403.2239, 2014.

1/--pages