1 Introduction
Blind superresolution of point sources is the problem of simultaneously estimating locations and amplitudes of point sources and point spread functions from lowresolution measurements. Such problem arises in various applications, including singlemolecule imaging
[13], medical imaging [12], multiuser communication system [10] and so on.Under certain subspace assumption and applying the lifting technique, blind superresolution can be cast as a matrix recovery problem. Recent works in [7, 17, 9, 15, 5] exploit the intrinsic structures of data matrix and propose different convex relaxation methods for such problem. Theoretical guarantees for these methods have been established. However, due to the limitations of convex relaxation, all of these methods do not scale well to the high dimensional setting.
In contrast to convex relaxation, a nonconvex recovery method is proposed in this paper based on the Vectorized Hankel Lift [5] framework. More precisely, harnessing lowrank structure of vectorized Hankel matrix corresponding to the signal in terms of the BurerMonteiro factorization, we develop a projected gradient descent algorithm, named PGDVHL, to directly recover the low rank factors. We show that such a simple algorithm possesses a remarkable reconstruction ability. Moreover, our algorithm started from an initial guess converges linearly to the target matrix under the similar sample complexity as convex approaches.
The rest of this paper is organized as follows. We begin with the problem formulation and describe the proposed algorithm in Section 2. Section 3 provides a convergence analysis of PGDVHL. Numerical experiments to illustrate the performance of PGDVHL are provided in Section 4. Section 5 concludes this paper.
2 Algorithm
2.1 Problem formulation
The point source signal can be represented as a superposition of spikes
where and are the amplitude and location of th point source respectively. Let be the unknown point spread functions. The observation is a convolution between and ,
After taking Fourier transform and sampling, we obtain the measurements as
(2.1) 
Let be a vector corresponding to the th unknown point spread function. The goal is to estimate as well as from (2.1). Since the number of unknowns is larger than , this problem is an illposed problem without any additional assumptions. Following the same route as that in [7, 17, 9, 5], we assume belong to a known subspace spanned by the columns of , i.e.,
Then under the subspace assumption and applying the lifting technique [1], the measurements (2.1) can be rewritten as a linear observations of :
(2.2) 
where is the th row of , is the th standard basis of , and is the vector defined as
The measurement model (2.2) can be rewritten succinctly as
(2.3) 
where is the linear operator. Therefore, blind superresolution can be cast as the problem of recovering the data matrix from its linear measurements (2.3).
Let be the vectorized Hankel lifting operator which maps a matrix into an matrix,
where is the th column of and . It is shown that is a rank matrix [5] and thus the matrix admits low rank structure when . It is natural to recover by solving the constrained least squares problem
To introduce our algorithm, we assume that is incoherent which is defined as below.
Assumption 2.1.
Let
be the singular value decomposition of
, where and . Denote , where is the th block of for . The matrix is incoherence if and obey thatfor some positive constant .
Remark 2.1.
Letting and , we have
where . Since that target data matrix is incoherence and the low rank structure of the vectorized Hankel matrix can be promoted by , it is natural to recover the low rank factors of the ground truth matrix by solving an optimization problem in the form of
(2.4) 
where , is the MoorePenrose pseudoinverse of obeying that , the second term in objective guarantees that admits vectorized Hankel structure. Last term penalizes the mismatch between and , which is widely used in rectangular low rank matrix recovery [16, 19, 6]. The convex feasible set is defined as follows
(2.5) 
where is the th block of , and be two absolute constants such that and .
2.2 Projected gradient descent method
Inspired by [2], we design a projected gradient descent method for the problem (2.1), which is summarized in Algorithm 1.
The initialization involves two steps: (1) computes the best rank approximation of via one step hard thresholding , where is the adjoint of ;(2) projects the low rank factors of best rank approximated matrix onto the set . Given a matrix , the projection onto , denoted by , has a closed form solution:
for and
for . Let be the current estimator. The algorithm updates along gradient descent direction with step size , followed by projection onto the set . The gradient of is computed with respect to Wirtinger calculus given by where
To obtain the computational cost of , we first introduce some notations. Let be the Hankel operator which maps a vector into an matrix,
where is the th entry of . The adjoint of , denoted by , is a linear mapping from to . It is known that the computational complexity of both and is flops [3]. Moreover, the authors in [5] show that , where is a matrix constructed by stacking all on top of one another, and is a permutation matrix. Therefore we can compute and by using flops. Thus the implementation of our algorithm is very efficient and the main computational complexity in each step is .
3 Main results
In this section, we provide an analysis of PGDVHL under a random subspace model.
Assumption 3.1.
The column vectors of are independently and identically drawn from a distribution which satisfies the following conditions
(3.1)  
(3.2) 
Remark 3.1.
Assumption 3.1 is a standard assumption in blind superresolution [7, 17, 9, 14, 5], and holds with by many common random ensembles, for instance, the components of
are Rademacher random variables taking the values
with equal probability or
is uniformly sampled from the rows of a Discrete Fourier Transform (DFT) matrix.Now we present the main result of the paper.
Theorem 3.1.
Remark 3.2.
Remark 3.3.
Theorem 3.1 implies that PGDVHL converges to with a linear rate. Therefore, after iterations, we have . Given the iterates returned by PGDVHL, we can estimate by .
Remark 3.4.
Once the data matrix is recovered, the locations can be computed from by MUSIC algorithm and the weights can be estimated by solving an overdetermined linear system [5].
4 Numerical simulations
In this section, we provide numerical results to illustrate the performance of PGDVHL. The locations of the point source signal is randomly generated from and the amplitudes are selected to be , where is uniformly sampled from and is uniformly sampled from . The coefficients are i.i.d. sampled from standard Gaussian with normalization. The columns of are uniformly sampled from the DFT matrix. The stepsize of PGDVHL is chosen via backtracking line search and PGDVHL will be terminated if is met or a maximum number of iterations is reached. We repeat 20 random trials and record the probability of successful recovery in our tests. A trial is declared to be successful if .
The first experiment studies the recovery ability of PGDVHL through the framework of phase transition and we compare it with two convex recovery methods: VHL
[5] and ANM [17]. Both VHL and ANM are solved by CVX [8]. The tests are conducted with and the varied and . Figure 1(a), 1(c) and 1(e) show that phase transitions of VHL, ANM and PGDVHL when the locations of point sources are randomly generated, and Figure 1(b), 1(d) and 1(f) illustrate the phase transitions of VHL, ANM and PGDVHL when the separation condition is imposed. It can be seen that PGDVHL is less sensitive to the separation condition than ANM and has a higher phase transition curve than VHL.In the second experiment, we study the phase transition of PGDVHL when one of and is ﬁxed. Figure 2(a) indicates a approximately linear relationship between and for the successful recovery when the number of point sources is ﬁxed to be . The same linear relationship between and can be observed when the dimension of the subspace is ﬁxed to be , see Figure 2(b). Therefore there exists a gap between our theory and empirical observation and we leave it as future work.
In the third simulation, we investigate the convergence rate of PGDVHL for different . The results are shown in Figure 3. The axis denotes and the axis represents the iteration number. It can be clearly seen that PGDVHL converges linearly as shown in our main theorem.
5 Conclusion
In this paper, we propose an efficient algorithm named PGDVHL towards recovering low rank matrix in blind superresolution. Our theoretical analysis shows that the proposed algorithm converges to the target matrix linearly. This is also demonstrated by our numerical simulations.
References
 [1] (2013) Blind deconvolution using convex programming. IEEE Transactions on Information Theory 60 (3), pp. 1711–1732. Cited by: §2.1.
 [2] (2018) Spectral compressed sensing via projected gradient descent. SIAM Journal on Optimization 28 (3), pp. 2625–2653. Cited by: §2.2.
 [3] (2019) Fast and provable algorithms for spectrally sparse signal reconstruction via lowrank Hankel matrix completion. Applied and Computational Harmonic Analysis 46 (1), pp. 94–121. Cited by: §2.2.
 [4] (2009) Exact matrix completion via convex optimization. Foundations of Computational Mathematics 9 (6), pp. 717. Cited by: Remark 2.1.
 [5] (2020) Vectorized hankel lift: a convex approach for blind superresolution of point sources. arXiv preprint arXiv:2008.05092. Cited by: §1, §1, §2.1, §2.1, §2.2, Remark 2.1, Remark 3.1, Remark 3.2, Remark 3.4, §4.
 [6] (2019) Nonconvex optimization meets lowrank matrix factorization: an overview. IEEE Transactions on Signal Processing 67 (20), pp. 5239–5269. Cited by: §2.1.
 [7] (2016) Guaranteed blind sparse spikes deconvolution via lifting and convex optimization. IEEE Journal of Selected Topics in Signal Processing 10 (4), pp. 782–794. Cited by: §1, §2.1, Remark 3.1.
 [8] (201403) CVX: matlab software for disciplined convex programming, version 2.1. http://cvxr.com/cvx. Cited by: §4.
 [9] (2019) Atomic norm denoising for complex exponentials with unknown waveform modulations. IEEE Transactions on Information Theory 66 (6), pp. 3893–3913. Cited by: §1, §2.1, Remark 3.1.
 [10] (2006) Lowcomplexity blind synchronization and demodulation for (ultra) wideband multiuser ad hoc access. IEEE Transactions on Wireless communications 5 (7), pp. 1930–1941. Cited by: §1.
 [11] (2021) Fast blind superresolution of point sources via projected gradient descent. Under Preparation. Cited by: Remark 3.2.
 [12] (2015) Accelerated NMR spectroscopy with lowrank reconstruction. Angewandte Chemie International Edition 54 (3), pp. 852–854. Cited by: §1.
 [13] (2012) Optimal 3D singlemolecule localization for superresolution microscopy with aberrations and engineered point spread functions. Proceedings of the National Academy of Sciences 109 (3), pp. 675–679. Cited by: §1.
 [14] (2018) Blind twodimensional superresolution and its performance guarantee. arXiv preprint arXiv:1811.02070. Cited by: Remark 3.1.
 [15] (2021) Mathematical theory of atomic norm denoising in blind twodimensional superresolution. IEEE Transactions on Signal Processing 69, pp. 1681–1696. Cited by: §1.

[16]
(2016)
Lowrank solutions of linear matrix equations via procrustes flow.
In
International Conference on Machine Learning
, pp. 964–973. Cited by: §2.1.  [17] (2016) Superresolution of complex exponentials from modulations with unknown waveforms. IEEE Transactions on Information Theory 62 (10), pp. 5809–5830. Cited by: §1, §2.1, Remark 3.1, §4.
 [18] (2018) Multichannel Hankel matrix completion through nonconvex optimization. IEEE Journal of Selected Topics in Signal Processing 12 (4), pp. 617–632. Cited by: Remark 2.1.
 [19] (2016) Convergence analysis for rectangular matrix completion using burermonteiro factorization and gradient descent. arXiv preprint arXiv:1605.07051. Cited by: §2.1.
Comments
There are no comments yet.