Compressive Sampling using Annihilating Filter-based Low-Rank Interpolation
1 While the recent theory of compressed sensing provides an opportunity to overcome the Nyquist limit in recovering sparse signals, a solution approach usually takes a form of inverse problem of the unknown signal, which is crucially dependent on specific signal representation. In this paper, we propose a drastically different two-step Fourier compressive sampling framework in continuous domain that can be implemented as a measurement domain interpolation, after which a signal reconstruction can be done using classical analytic reconstruction methods. The main idea is originated from the fundamental duality between the sparsity in the primary space and the low-rankness of a structured matrix in the spectral domain, which shows that a low-rank interpolator in the spectral domain can enjoy all the benefit of sparse recovery with performance guarantees. Most notably, the proposed low-rank interpolation approach can be regarded as a generalization of recent spectral compressed sensing to recover large class of finite rate of innovations (FRI) signals at near optimal sampling rate. Moreover, for the case of cardinal representation, we can show that the proposed low-rank interpolation will benefit from inherent regularization and the optimal incoherence parameter. Using the powerful dual certificates and golfing scheme, we show that the new framework still achieves the near-optimal sampling rate for general class of FRI signal recovery, and the sampling rate can be further reduced for the class of cardinal splines. Numerical results using various type of FRI signals confirmed that the proposed low-rank interpolation approach has significant better phase transition than the conventional CS approaches.
Jong Chul Ye, Ph.D
Department of Bio and Brain Engineering
Korea Adv. Inst. of Science and Technology (KAIST)
373-1 Guseong-Dong, Yuseong-Gu, Daejon 305-701, Korea
Compressed sensing or compressive sampling (CS) theory [1, 2, 3] addresses the accurate recovery of unknown sparse signals from underdetermined linear measurements. In particular, a Fourier CS problem, which recovers unknown signals from sub-sampled Fourier measurements, has many important applications in imaging applications such as magnetic resonance imaging (MRI), X-ray computed tomography (CT), optics, and so on. Moreover, this problem is closely related to the classical harmonic retrieval problem that computes the amplitudes and frequencies at off the grid locations of a superposition of complex sinusoids from their consecutive or bunched Fourier samples. Harmonic retrieval can be solved by various methods including Prony’s method, and matrix pencil algorithm . These methods were proven to succeed at the minimal sample rate in the noiseless case, because it satisfies an algebraic condition called the full spark (or full Kruskal rank) condition  that guarantees the unique identification of the unknown signal. Typically, when operating at the critical sample rate, these method are not robust to perturbations in the measurements due to the large condition number.
Accordingly, to facilitate robust reconstruction of off the grid spectral components, CS algorithms from non-consecutively sub-sampled Fourier measurements are required. The scheme is called spectral compressed sensing, which is also known as compressed sensing off the grid, when the underlying signal is composed of Diracs. Indeed, this has been developed with a close link to the recent super-resolution researches [7, 8, 9, 10]. For example, Candes and/or Fernandez-Granda [8, 11] showed that if the minimum distance of the Diracs is bigger than where denotes the cut-off frequency of the measured spectrum, then a simple convex optimization can solve the locations of the Diracs. Under the same minimum distance condition, Tang [9, 12] proposed an atomic norm minimization approach for the recovery of Diracs from random spatial and Fourier samples, respectively. Unlike these direct signal recovery methods, Chen and Chi  proposed a two-step approach consisting of interpolation followed by a matrix pencil algorithms. In addition, they provided performance guarantees at near optimal sample complexity (up to a logarithmic factor). One of the main limitations of these spectral compressed sensing approaches is, however, that the unknown signal is restricted to a stream of Diracs. The approach by Chen and Chi  is indeed a special case of the proposed approach, but they did not realize its potential for recovering much wider class of signals.
Note that the stream of Diracs is a special instance of a signal model called the signals with the finite rate of innovation (FRI) [13, 14, 15]. Originally proposed by Vetterli et al , the class of FRI signals includes a stream of Diracs, a stream of differentiated Diracs, non-uniform splines, piecewise smooth polynomials, and so on. Vetterli et al [13, 14, 15] proposed time-domain sampling schemes of these FRI signals that operate at the rate of innovation with a provable algebraic guarantee in the noise-free scenario. Their reconstruction scheme estimates an annihilating filter that cancels the Fourier series coefficients of a FRI signal at consecutive low-frequencies. However, due to the time domain data acquisition, the equivalent Fourier domain measurements are restricted to a bunched sampling pattern similar to the classical harmonic retrieval problems.
Therefore, one of the main aims of this paper is to generalize the scheme by Verterli et al [13, 14, 15] to address Fourier CS problems that recover general class of FRI signals from irregularly subsampled Fourier measurements. Notably, we prove that the only required change is an additional Fourier domain interpolation step that estimates missing Fourier measurements. More specifically, for general class FRI signals introduced in [13, 14, 15], we show that there always exists a low-rank Hankel structured matrix associated with the corresponding annihilating filter. Accordingly, their missing spectral elements can be interpolated using a low-rank Hankel matrix completion algorithm. Once a set of Fourier measurements at consecutive frequencies are interpolated, a FRI signal can be reconstructed using conventional methods including Prony’s method and matrix pencil algorithms as done in [13, 14, 15]. Most notably, we show that the proposed Fourier CS of FRI signals operates at a near optimal rate (up to a logarithmic factor) with provable performance guarantee. Additionally, thanks to the inherent redundancies introduced by CS sampling scheme, the subsequent step of retrieving a FRI signal becomes much more stable.
While a similar low-rank Hankel matrix completion approach was used by Chen and Chi , there are several important differences. First, the low-rankness of the Hankel matrix in  was shown based on the standard Vandermonde decomposition, which is true only when the underlying FRI signal is a stream of Diracs. Accordingly, in case of differentiated Diracs, the theoretical tools in  cannot be used. Second, when the underlying signal can be converted to a stream of Diracs or differentiated Diracs by applying a linear transform that acts as a diagonal operator (i.e., element-wise multiplication) in the Fourier domain, we can still construct a low rank Hankel matrix from the weighted Fourier measurements, whose weights are determined by the spectrum of the linear operator. For example, a total variation (TV)-sparse signal is a stream of Diracs after the differentiation, and piecewise smooth polynomials becomes a stream of differentiated Diracs by applying a differential operator. Finally, the advantage of the proposed approach becomes more evident when we model the unknown signal using cardinal L-splines . In cardinal L-splines, the discontinuities occur only on an integer grid, which is a reasonable model to acquire signals of high but finite resolution. Then, we can show that the discretization using cardinal splines makes the reconstruction significantly more stable in the existence of noise to measurements, and the logarithmic factor as well as the incoherence parameter for the performance guarantees can be further improved.
It is important to note that the proposed low-rank interpolation approach is different from classical compressed sensing approaches which regard a sampling problem as an inverse problem and whose focus is to directly recover the unknown signal. Rather, the proposed approach is more closely related to the classical sampling theory, where signal sampling step is decoupled from a signal recovery algorithm. For example, in the sampling theory for signals in the shift-invariant spaces [17, 18], the nature of the signal sampling can be fully taken care of as a digital correction filter, after which signal recovery is performed by convolution with a reconstruction filter (see Fig. 1(a)). Similarly, by introducing a low rank interpolator, the proposed scheme in Fig. 1(b) fully decouples the signal recovery from sampling step as a separate layer that can be optimized independently. This is because the same low-rank interpolator will successfully complete missing measurements, regardless of whether the unknown signal is either a stream of Diracs or a stream of differentiated Diracs. In the subsquent step, analytic reconstruction methods such as Prony’s method and matrix pencil algorithms can identify the signal model as done in [13, 14, 15].
The proposed two-layer approach composed of Fourier domain interpolation and the analytic reconstruction is very useful in real-world applications, because the low-rank interpolator can be added as a digital correction filter to existing systems, where the second step is already implemented. Moreover, in many biomedical imaging problems such as magnetic resonance imaging (MRI) or X-ray computed tomography (CT), an accurate interpolation to fully sampled Fourier data gives an important advantage of utilising fully established mathematical theory of analytic reconstruction. In addition, classical preprocessing techniques for artifact removal have been developed by assuming fully sampled measurements, so these steps can be readily combined with the proposed low-rank interpolation approaches. The superior advantages of the proposed scheme have been demonstrated in various biomedical imaging and image processing applications such as compressed sensing MRI [19, 20], MR artifact correction[21, 22], image inpainting , super-resolution microscopy , image denoising , and so on, which clearly confirm the practicality of the new theory.
Nonetheless, it is remarkable that the proposed two-layer approach using a low-rank interpolation achieves near optimal sample rate while universally applying to different signal models of the same order (e.g., stream of Diracs and stream of differentiated Diracs). Moreover, it may look mysterious that no explicit form of the minimum separation distance as required in Fernandez-Granda [8, 11] and Tang [9, 12] is not shown in the performance guarantee. However, the proposed method is not free of limitations. Specifically, we will show that the incoherence parameter in our performance guarantees is dependent upon the type of unknown signals as well as the minimum separation between the successive spikes. The similarity and differences of our results from the existing theory [8, 11, 9, 12] and the origin of the differences will be also discussed.
This paper consists of followings. Section II first discusses the main results that relate an annihilating filter and a low-rank Hankel structured matrix, and provides the performance guarantees of low-rank structured matrix completion, which will be used throughout the paper. Section III then discusses the proposed low-rank interpolation theory for recovery of FRI signals, which is followed by the low rank interpolation for the case of cardinal L-splines in Section IV. Section V explains algorithmic implementation. Numerical results are then provided in Section VI, which is followed by conclusion in Section VII, respectively.
Ii Main Results
A Hankel structured matrix generated from an -dimensional vector has the following structure:
where is called a matrix pencil parameter. We denote the space of this type of Hankel structure matrices as .
An wrap-around Hankel matrix generated from an -dimensional vector is defined as:
Note that wrap-around Hankel matrix can be considered as a Hankel matrix of -element augumented vector from with the periodic boundary expansion:
We denote the space of this type of wrap-around Hankel structure matrices as .
Ii-B Annihilating Filter-based Low-Rank Hankel Matrix
The Fourier CS problem of our interest is to recover the unknown signal from the Fourier measurement:
Without loss of generality, we assume that the support of is . Then, the sampled Fourier data at the Nyquist rate is defined by
We also define a length -annihilating filter for that satisfies
Suppose that the filter is the minimum length annihilating filter. Then, for any tap filter , it is easy to see that the following filter with taps is also an annihilating filter for :
because . The matrix representation of (5) is given by
where denotes a vector that reverses the order of the elements in
Accordingly, by choosing such that and defining an -dimensional vector composed of sampled Fourier data at the Nyquist rate as:
we can construct the following matrix equation:
where the Hankel structure matrix is constructed as
Then, we can show the following key result:
Let denote the minimum size of annihilating filters that annihilates sampled Fourier data . Assume that . Then, for a given Hankel structured matrix constructed in (10), we have
where denotes a matrix rank.
See Appendix -B. ∎
Ii-C Performance Guarantees for Structured Matrix Completion
Let be a multi-set consisting of random indices from such that . While the standard CS approaches directly estimate from , here we propose a two-step approach by exploiting Theorem II.1. More specifically, we first interpolate for all from the sparse Fourier samples, and the second step then applies the existing spectral estimation methods to estimate as done in [13, 14, 15]. Thanks to the low-rankness of the associated Hankel matrix, the first step can be implemented using the following low-rank matrix completion:
where is the projection operator on the sampling location . Therefore, the remaining question is to verify whether the low-rank matrix completion approach (12) does not compromise any optimality compared to the standard Fourier CS, which is the main topic in this section.
The low-rank matrix completion problem in (12) is non-convex, which is difficult to analyze. Therefore, to provide a performance guarantee, we resort to its convex relaxation using the nuclear norm. Chen and Chi  provided the first performance guarantee for robust spectral compressed sensing via structured matrix completion by nuclear norm minimization and extended the result to general low-rank Hankel/Toeplitz matrix completion [10, Theorem 4]. However, parts of their proof (e.g., [10, Appendix H]) critically depend on the special structure given in the standard Vandermonde decomposition. Furthermore, they also used the standard incoherence condition, which is neither assumed nor applied by their incoherence condition [10, eq. (27)]. Therefore, unlike their claim, the main results in , in its current forms, apply only to the spectral compressed sensing. Here, we elaborate on their results so that the performance guarantees apply to the general structured low-rank matrix completion problems which will be described in subsequent sections.
Recall that the notion of the incoherence plays a crucial role in matrix completion and structured matrix completion. We recall the definitions using our notations. Suppose that is a rank- matrix whose SVD is with and , respectively. is said to satisfy the standard incoherence condition with parameter if
where denotes the appropriate size standard coordinate vector with 1 on the -th elements and zeros elsewhere.
To deal with two types of Hankel matrices simultaneously, we define a linear lifting operator that lifts a vector to a structured matrix in a higher dimensional space. For example, the dimension of is given as and for a lifting to a Hankel matrix in , whereas and for a lifting to a wrap-around Hankel matrix . A linear lifting operator is regarded as a synthesis operator with respect to basis :
where the specific form of the basis for the case of for and will be explained in Appendix -C.
Then, the completion of a low rank structured matrix from the observation of its partial entries can be done by minimizing the nuclear norm under the measurement fidelity constraint as follows:
where denotes the matrix nuclear norm. Then, we have the following main result.
Let be a multi-set consisting of random indices where ’s are i.i.d. following the uniform distribution on . Suppose correspond to one of the structured matrices in and . Suppose, furthermore, that is of rank- and satisfies the standard incoherence condition in (13) with parameter . Then there exists an absolute constant such that is the unique minimizer to (14) with probability , provided
where if each images of has the wrap-around property; , otherwise, and .
See Appendix -E. ∎
Note that Theorem II.2 provides a generalized version of performance guarantee compared to the previous work . Specifically, Theorem II.2 holds for other structured matrices, if the associated basis matrix satisfies the specific condition described in detail in Eq. (A.17). In addition, for the case of signals with wrap-around Hankel matrix (which will be explained later), the log exponentional factor in (15) becomes 2, which reduces the sampling rate. This sampling rate reduction is novel and was not observed in the previous work .
Next, we consider the recovery of from its partial entries with noise. Let denote a corrupted version of . The unknown structured matrix can be estimated via
Then, we have the following stability guarantee:
Proof of Theorem ii.3.
Iii Guaranteed Reconstruction of FRI Signals
The explicit derivation of the minimum length finite length annihilating filter was one of the most important contributions of the sampling theory of FRI signals [13, 14, 15]. Therefore, by combing the results in the previous section, we can provide performance guarantees for the recovery of FRI signals from partial Fourier measurements.
Iii-a Spectral Compressed Sensing: Recovery of Stream of Diracs
Consider the periodic stream of Diracs described by the superposition of impulses
Then, the discrete Fourier data are given by
One of the important contributions of this section is to show that the spectral compressed sensing can be equivalently explained using annihilating filter-based low-rank Hankel matrix. Specifically, for the stream of Diracs, the minimum length annihilating filter has the following z-transform representation :
For a given stream of Diracs in Eq. (17), denotes the noiseless discrete Fourier data in (8). Suppose, furthermore, is given by . Let is a multi-set consisting of random indices where ’s are i.i.d. following the uniform distribution on . Then, there exists an absolute constant such that is the unique minimizer to (III-A) with probability , provided
This result appears identical to that of Chen and Chi . However, they explicitly utilized the standard Vandermonde decomposition. On the contrary, the annihilating filter-based construction of low-rank Hankel matrix is more general that can cover all FRI signal models as will be shown later.
For the noisy measurement, we interpolate the missing Fourier data using the following low-rank matrix completion:
Iii-B Stream of Differentiated Diracs
Another important class of FRI signal is a stream of differentiated Diracs:
where denotes the -th derivative of Diracs in the distributions sense. Thus, its Fourier transform is given by
whose discrete samples are given by
Then, there exists an associated minimum length annihilating filter whose z-transform is given by:
where . Therefore, we can provide the following performance guarantees:
For a given stream of differentiated Diracs in Eq. (24), denotes the noiseless discrete Fourier data in (8). Suppose, furthermore, is given by . Let be a multi-set consisting of random indices where ’s are i.i.d. following the uniform distribution on . Then, there exists an absolute constant such that is the unique minimizer to (III-A) with probability , provided
Iii-C Non-uniform Splines
Note that signals may not be sparse in the image domain, but can be sparsified in a transform domain. Our goal is to find a generalized framework, whose sampling rate can be reduced down to the transform domain sparsity level. Specifically, the signal of our interest is a non-uniform spline that can be represented by :
and is a continuous sparse innovation:
For example, if the underlying signal is piecewise constant, we can set as the first differentiation. In this case, corresponds to the total variation signal model. Then, by taking the Fourier transform of (29), we have
Accordingly, the same filter whose z-transform is given by (19) can annihilate the discrete samples of the weighted spectrum , and the Hankel matrix from the weighted spectrum satisfies the following rank condition:
Thanks to the low-rankness, the missing Fourier data can be interpolated using the following matrix completion problem:
or, for noisy Fourier measurements ,
where denotes the Hadamard product, and and denotes the vectors composed of full samples of and , respectively. After solving , the missing spectral data can be obtained by dividing by the weight, i.e. assuming that , where is the estimated Fourier data using . As for the sample at the spectral null of the filter , the corresponding elements should be included as measurements.
Now, we can provide the following performance guarantee:
For a given non-uniform splines in Eq. (29), denotes the noiseless discrete Fourier data in (8). Suppose, furthermore, is given by and be a multi-set consisting of random indices where ’s are i.i.d. following the uniform distribution on . Then, there exists an absolute constant such that is the unique minimizer to (III-C) with probability , provided
Iii-D Piecewise Polynomials
A signal is a periodic piecewise polynomial with pieces each of maximum degree if and only if its derivative is a stream of differentiated Diracs given by
In this case, the corresponding Fourier transform relationship is given by
whose filter length is given by . Therefore, we can provide the following performance guarantee:
For a given piecewise smooth polynomial in Eq. (37), let denotes the discrete spectral samples of with . Suppose, furthermore, is given by and be a multi-set consisting of random indices where ’s are i.i.d. following the uniform distribution on . Then, there exists an absolute constant such that is the unique minimizer to (III-C) with probability , provided
Iii-E Incoherence and the Minimum Separation
Note that the proposed low-rank interpolation achieves near optimal sample rate while universally applying to different FRI signal models of the same order (e.g., stream of Diracs and stream of differentiated Diracs). Moreover, even though the concept of the minimum separation between the successive spikes was essential for the performance guarantee of super-resolution in Candes and Fernandez-Granda [8, 11], Tang [9, 12], etc, similar expression is not observable in Theorem III.4-Theorem III.11. This looks mysterious. Therefore, the main goal of this section is to show that these information are still required but hidden in the incoherence parameter .
Note that the proof in Theorem II.1 implies that the explicit form given by
is a necessary and sufficient condition to have the low-rank Hankel matrix. Here, for FRI signals. Furthermore, for a Hankel matrix constructed using the signal model in (41), there exist an exponentional decomposition of Hankel matrix using confluent Vandermonde matrix . Specifically, define a confluent Vandermonde matrix (resp. ):
where the element of the sub-matrix is given by
where and are the confluent Vandermonde matrices and is a block diagonal matrix given by
Accordingly, we can derive the following upper bound of the standard coherence:
Note that this is an extension of the approach in [10, Appendix C, Section III.A] which tried to bound the mutual coherence for the standard Vandermonde decomposition, (i.e. ) by a small number. Specifically, for the cases of random frequency locations or small perturbation off the grid, they showed that the incoherence parameters become small . However, the dependency of on the minimum separation was not explicit in their discussion.
Recently, Moitra  discovered a very intuitive relationship between the least/largest singular values of Vandermonde matrix and the minimum separation distance Specifically, by making novel connections between extremal functions and the spectral properties of Vandermonde matrices , Moitra showed that if , then the least singular value is bounded as
If applied in our problem involving and , the resulting upper bound of the coherence parameter for standard Vandermonde matrix is given by
which approaches to one with sufficiently large . Eq. (47) is obtained because the matrix pencil size gives the optimal trade-off between and , and owing to . Note that this coincides with the minimum separation in Fernandez-Granda [8, 11] and Tang [9, 12]. However, compared to these approaches [8, 11, 9, 12] that require the minimum separation as a hard constraint, our approach requires it as a soft oversampling factor in terms of the incoherence parameter .
Then, where is the difference originated ? We argue that this comes from different uses of interpolation functions. Specifically, in Candes, Fernandez-Granda [8, 11] and Tang [9, 12], dual polynomial function that interpolates the sign at the singularity locations should be found to construct a dual certificate. On other hand, in the proposed annihilating filter based approach, the interpolating function is a smooth function that has zero-crossings at the singularity locations. To see this, let is an annihilating filter that annihilates . Then there exists an annihilating function such that
so whenever . The construction of the annihilating function is extremely easy and can be readily obtained by the multiplications of sinusoids (for example, to null out -periodic stream of Diracs within , we set ). Moreover, this approach can be easily extended to have multiple roots, which is required for differentiated Diracs. We believe that the “soft constraint” originated from annihilating function is one of the key ingredients that enables recovery of general FRI signals which was not possible by the existing super-resolution methods [8, 11, 9, 12].
The derivation of the least singular value for the confluence Vandermonde matrix have been also an important topic of researches [33, 31, 30, 34, 35, 36, 37]. In general, it will also depend on the minimum separation distance . However, the explicit tight bound of the least singular value is not available in general, so we leave this for future work.
Iii-F Recovery of Continuous Domain FRI Signals After Interpolation
Regardless of the unknown signal type (the stream of Diracs or a stream of differentiated Diracs), note that an identical low-rank interpolator can be used. Once the spectrum is fully interpolated, in the subsequent step, Prony’s method and matrix pencil algorithm can identify the signal model from the roots of the estimated annihilator filter as done in [13, 14, 15]. Accordingly, our robustness guarantees on the low-rank matrix entries can be translated in terms of the actual signal that is recovered (for example, on the support or amplitudes of the spike in the case of recovery of spike superpositions). In fact, this has been also an active area of researches [33, 31, 30, 32, 38, 29], and we again exploit these findings for our second step of signal recovery. For example, see Moitra  for more details on the error bound for the case of modified matrix pencil approach for recovery of Diracs. In addition, Batenkov has recently generalized this for the recovery of general signals in (41) . The common findings are that the estimation error for the location parameter and the magnitude are bounded by the condition number of the confluent Vandermonde matrix as well as the minimum separation distance . Moreover, matrix pencil approaches such as Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) method  is shown to stably recovery the locations [38, 33].
Here, we briefly review the matrix pencil approach for the signal recovery [29, 40]. Specifically, for a given confluent Vandermonde matrix , let be the matrix extracted from by deleting the last row. Similarly, let be the matrix extracted from by deleting the first row. Then, and span the same signal subspace and
where is the block diagonal matrix