close

Вход

Log in using OpenID

embedDownload
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
Пожаловаться на содержимое документа