Structured low-rank recovery of piecewise constant signals with performance guarantees
We derive theoretical guarantees for the exact recovery of piecewise constant two-dimensional images from a minimal number of non-uniform Fourier samples using a convex matrix completion algorithm. We assume the discontinuities of the image are localized to the zero level-set of a bandlimited function, which induces certain linear dependencies in Fourier domain, such that a multifold Toeplitz matrix built from the Fourier data is known to be low-rank. The recovery algorithm arranges the known Fourier samples into the structured matrix then attempts recovery of the missing Fourier data by minimizing the nuclear norm subject to structure and data constraints. This work adapts results by Chen and Chi on the recovery of isolated Diracs via nuclear norm minimization of a similar multifold Hankel structure. We show that exact recovery is possible with high probability when the bandlimited function describing the edge set satisfies an incoherency property. Finally, we demonstrate the algorithm on the recovery of undersampled MRI data.
Structured low-rank recovery of piecewise constant signals with performance guarantees
|Greg Ongie, Sampurna Biswas, and Mathews Jacob††thanks: This work is supported by grants NSF CCF-1116067, ACS RSG-11-267-01-CCE, and ONR-N000141310202.|
|Department of Mathematics, University of Iowa, IA, USA.|
|Department of Electrical and Computer Engineering, University of Iowa, IA, USA.|
Index Terms— Structured Low-Rank Matrix Completion, Annihilating Filter Method, Finite Rate of Innovation, Compressed Sensing, MRI.
The recovery of a linear combination of exponentials from their few uniform samples is a classical problem in signal processing with extensive applications. Prony’s method, or one of its robust variants, attempts to recover the signal by estimating an annihilating polynomial whose zeros correspond to the frequency of the exponentials. The finite rate of innovation (FRI) framework  extended these methods to recover more general signals that reduce to a sparse linear combination of Dirac delta functions under an appropriate transformation (e.g. differential operators, convolution). Recently, several authors have further extended FRI methods to recover such signals from their non-uniform Fourier samples [2, 3, 4, 5, 6] by exploiting the low-rank structure of an enhanced matrix (e.g. Hankel matrix in 1-D). Performance guarantees do exist when the transform is an identity and when the Diracs are well-separated .
The above signal models have limited flexibility in exploiting the extensive additional structure present in many multidimensional imaging problems. Specifically, the edges in multidimensional images are connected and can be modeled as smooth curves or surfaces. We have recently introduced a novel framework to recover piecewise polynomial images, whose edges are localized to smooth curves, from their uniform [7, 8] and non-uniform  Fourier samples; this work generalizes a recent extension of the FRI framework to curves . We model the piecewise smooth signal as having partial derivatives that vanish outside the zero level-set of a bandlimited function. This relation translates to an annihilation condition involving the uniform Fourier samples of the partial derivatives, which can be compactly represented as the multiplication of a specific structured matrix with the Fourier coefficients of the bandlimited function. Our earlier work has shown that the structured matrix is low-rank, and we used this property to recover the signal from its non-uniform Fourier samples with good performance. Efficient algorithms that work on the original signal samples rather than the structured high-dimensional matrix also were introduced . We observe the signal models in [2, 3, 5] do not include the class of signals considered in this work.
The main focus of this work is to introduce theoretical guarantees on the recovery of piecewise constant signals, whose discontinuities are localized to zero level-sets of bandlimited functions, from non-uniform Fourier samples. Since such signals cannot be expressed as a finite linear combination of isolated Diracs, the recovery guarantees in  cannot be directly extended to our setting. Specifically, the theory in  relies heavily on a explicit factorization of the enhanced matrix (e.g Vandermonde factorization of a Hankel matrix in the 1-D case), which is only available when the number of discontinuities are finite and well separated. Instead, we give a new description of the row and column subspace of the structured matrix, which allow us to derive incoherence measures based solely on properties of the bandlimited function describing the edge set of the image.
2.1 Signal Model: 2-D Piecewise Constant Images
We consider the recovery of a piecewise constant function
where , and denotes the characteristic function of the set . We assume the Fourier samples specified by
are available at a subset of non-uniform locations belonging to a rectangular set of uniform sampling locations in Fourier domain.
We further assume that the edge-set of the image, specified by , to be the zero-set of a 2-D bandlimited trigonometric polynomial:
where the coefficients , and is a rectangular subset of . Here we assume is the unqiue minimal degree trigonometric polynomial such that , where the degree is defined by the rectangular dimensions of the Fourier support . We have shown in  that when is any bandlimited function that vanishes on the edgeset , the gradient satisfies the property
in the distributional sense. See Fig. 1 for an illustration when the number of regions . Note that among all functions bandlimited to in Fourier domain, is the only one in this class that satisfies (4). However, if we consider that are bandlimited to a larger rectangular set with , then we have shown that all satisfying (4) are a multiple of . The spatial domain relation (4) translates directly to the following vector annihilation relation in the Fourier domain:
Here for , and is any rectangular set on which the convolutions between and is valid. Note that the 2-D convolution between two filters support limited to rectangular sets and is support limited to the dilation of by , which we denote by . Since is bandlimited to , when using samples of within , we require ; see Fig. 2.
The Fourier domain annihilation relations (5) can be compactly represented in matrix form as
where , , are matrices corresponding to the discrete 2-D convolution of , , (omitting the irrelevant factor ) with a filter supported on , with output restricted to the index set . Here we use to denote the vectorized version of a filter . By our previous observations, the solutions to (6) are given by the Fourier coefficients of a multiple of the minimal polynomial. Hence if the filter support is larger than the minimal filter support , has a large nullspace and is low-rank. Specifically, in  we proved the following:
Suppose is built with filter size satisfying , then
where is the number of indices in and is the number of integer shifts of contained in .
2.2 Recovery from non-uniform Fourier samples
Since the matrix is low-rank, we propose to recover the signal from its noiseless Fourier samples as the convex optimization problem:
where denotes the nuclear norm, i.e. the sum of the singular values. To aid in our analysis, we will now reformulate the recovery of as a matrix recovery problem using projection operators in the lifted matrix domain. We define basis matrices , for all , where
for . Here is the set of locations of the matrix containing copies of . Note that the set forms an orthonormal basis for the space of structured matrices defined by the lifting . Specifically, for any set of coefficients we can expand the matrix as , where . Using these basis matrices, we define the following operators in the lifted domain:
Several authors have shown that the performance of low-rank matrix recovery by nuclear norm minimization is dependent on the incoherence of sampling basis with respect to the matrix to be to be recovered [11, 2]. Towards this end, we introduce two incoherence measures associated with defined solely in terms of the edge-set polynomial . In the following, we set to be the 2-D Dirichlet kernel supported on , i.e. the function such that for all and zero otherwise. For any collection of points , we define the Gram matrix by .
Let be a trigonometric polynomial bandlimited to (see (3)), and set . Define the incoherence measure by
where is the minimum singular value of .
Put in words, among all possible arrangements of points along the edge-set , we seek the arrangement that gives the best conditioning of the matrix , and call the resulting condition number . Intuitively, the optimal arrangement will maximize the minimum separation distance among the points, and can be thought of as a measure of this geometric property. In particular, having any edges that enclose a small area will give a high .
Additionally, our results rely on another incoherence measure related to properties of the gradient of :.
Let be a trigonometric polynomial bandlimited to (see (3)). Normalize such that . Define the incoherence measure by
where denotes the space of all trigonometric polynomials bandlimited to , and .
Note that will be large when has several zeros, or equivalently, when has several critical points. Since must contain a critical point in every region defined by the complement of the edge-set, will be large when the image has several distinct regions.
Now we state our main result:
Let be specified by (1), whose edge-set is described by the zero-set of the trigonometric polynomial bandlimited to (see (3)) with associated incoherence measures and . Let be an index set drawn uniformly at random within . Then there exists a universal constant such that the solution to (12) is exact () with probability exceeding , provided
where and .
Following the approach in , we can prove this result by constructing an approximate dual certificate using the well-known golfing scheme of . The adaptation of the proof in  to the measurement operator (11) is straightforward, and these details are omitted for brevity. The essential difference of the above result and  is the characterization of the incoherency measures. The approach in  relies on an explicit low-rank factorization of the lifted matrix in terms of Vandermonde-like matricies, which is not available in our setting. Instead, we characterize the row and column spaces of the matrix and use it to prove the above result.
2.3 Row and column spaces of and incoherence
Define and to be the orthogonal projections onto the column space and row space of , respectively, i.e. if is the rank- singular value decomposition then , . One can show it is possible to construct an approximate dual certificate with high probability [2, 11], provided we can uniformly bound the norms of the projections and . The following proposition shows these norms can be controlled by the incoherence measures introduced in (13) and (14)
Consider of rank corresponding to a piecewise constant function whose edge set coincides with the zero set of , let and be the incoherency measures of , and set . Then we have
The proof relies on the following basis representations for the row and column spaces of :
Choose any points , and define to be the collection of filters where each are the Fourier coefficients of the translated Dirichlet kernel . Then there exists a subset elements from that is a basis for the row space of .
Let be any basis of the row space of , and set . Then is a basis of the column space of .
In Fig. 3 we demonstrate the recovery of a synthetic piecewise constant phantom  ( sampling grid, single channel) from 20-fold variable density random Fourier samples using the structured low-rank matrix completion approach (7). We solve (7) using a singular value thresholding approach proposed in . The filter size was set to . Compared with total variation (TV) minimization, the proposed structured low-rank approach more accurately recovers the original piecewise constant regions.
Additionally, in Fig. 4 we demonstrate the structured low-rank approach on the recovery real MR data ( sampling grid with 4 receiver coils, coil-compressed to a single channel) from 2-fold uniform random undersampling, using a filter size of . Due to the problem size, the formulation (7) is difficult to solve via singular value thresholding. Instead, we make use of the recently proposed GIRAF algorithm  which solves an approximated version of (7). The result shows similar benefit over a TV regularized recovery in its ability to preserve fine details and strong edges.
4 Discussion and Conclusion
We derived a performance guarantee for the recovery of piecewise constant images from non-uniform Fourier samples by a structured matrix completion. This was achieved by adapting results in  to the case of a low-rank multifold Toeplitz structure with an additional weighting scheme. We also define new incoherence measures that rely only on properties of the minimal annihilating polynomial whose zero-set encodes the edges of the image.
While in the present work we only consider noiseless ideal samples, in future work we intend to derive guarantees for robust recovery in the presence of noise and model-mismatch. Additionally, it would be interesting to adapt our results to a wider variety of sampling distributions, and to identify the optimal sampling strategy for signals belonging to our image model.
-  M. Vetterli, P. Marziliano, and T. Blu, “Sampling signals with finite rate of innovation,” Signal Processing, IEEE Transactions on, vol. 50, no. 6, pp. 1417–1428, 2002.
-  Y. Chen and Y. Chi, “Robust spectral compressed sensing via structured matrix completion,” Information Theory, IEEE Transactions on, vol. 60, no. 10, pp. 6576–6601, 2014.
-  K. H. Jin, D. Lee, and J. C. Ye, “A general framework for compressed sensing and parallel MRI using annihilating filter based low-rank hankel matrix,” arXiv preprint arXiv:1504.00532, 2015.
-  ——, “A novel k-space annihilating filter method for unification between compressed sensing and parallel mri,” in International Symposium on Biomedical Imaging (ISBI), April 2015.
-  J. P. Haldar, “Low-rank modeling of local k-space neighborhoods (LORAKS) for constrained MRI.” Medical Imaging, IEEE Transactions on, vol. 33, no. 3, pp. 668–681, 2014.
-  G. Ongie and M. Jacob, “Recovery of piecewise smooth images from few fourier samples,” in Sampling Theory and Applications (SampTA), May 2015, pp. 543–547.
-  ——, “Super-resolution MRI using finite rate of innovation curves,” IEEE International Symposium on Biomedical Imaging (ISBI), April 2015.
-  ——, “Off-the-grid recovery of piecewise constant images from few fourier samples,” in arXiv:1510.00384, May 2015, pp. 543–547.
-  H. Pan, T. Blu, and P. L. Dragotti, “Sampling curves with finite rate of innovation,” Signal Processing, IEEE Transactions on, vol. 62, no. 2, 2014.
-  G. Ongie and M. Jacob, “A fast algorithm for structured low-rank matrix recovery with applications to undersampled MRI reconstruction,” IEEE International Symposium on Biomedical Imaging (ISBI), April 2016.
-  D. Gross, “Recovering low-rank matrices from few coefficients in any basis,” Information Theory, IEEE Transactions on, vol. 57, no. 3, pp. 1548–1566, 2011.
-  M. Guerquin-Kern, L. Lejeune, K. P. Pruessmann, and M. Unser, “Realistic analytical phantoms for parallel magnetic resonance imaging,” Medical Imaging, IEEE Transactions on, vol. 31, no. 3, pp. 626–636, 2012.