Abstract
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.
Keywords: Structured Low-Rank Matrix Completion, Annihilating Filter Method, Finite Rate of Innovation, Compressed Sensing, MRI
1. INTRODUCTION
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 [1] 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–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 [2].
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 [6] Fourier samples; this work generalizes a recent extension of the FRI framework to curves [9]. 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 [10]. 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 [2] cannot be directly extended to our setting. Specifically, the theory in [2] 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. THEORY
2.1. Signal Model: 2-D Piecewise Constant Images
We consider the recovery of a piecewise constant function
(1) |
where , and denotes the characteristic function of the set Ω. We assume the Fourier samples specified by
(2) |
are available at a subset of non-uniform locations k ∈ Θ belonging to a rectangular set of uniform sampling locations in Fourier domain.
We further assume that the edge-set of the image, specified by E := ∪i∂Ωi, to be the zero-set of a 2-D bandlimited trigonometric polynomial:
(3) |
where the coefficients , and Λ0 is a rectangular subset of Γ. Here we assume μ0 is the unqiue minimal degree trigonometric polynomial such that E = {μ0 = 0}, where the degree is defined by the rectangular dimensions of the Fourier support Λ0. We have shown in [8] that when μ(r) is any bandlimited function that vanishes on the edgeset E, the gradient ∇f = (∂xf, ∂yf) satisfies the property
(4) |
in the distributional sense. See Fig. 1 for an illustration when the number of regions N = 1. Note that among all functions bandlimited to Λ0 in Fourier domain, μ = μ0 is the only one in this class that satisfies (4). However, if we consider μ that are bandlimited to a larger rectangular set Λ1 with Λ0 ⊂ Λ1 ⊂ Γ, then we have shown that all μ satisfying (4) are a multiple of μ0 [8]. The spatial domain relation (4) translates directly to the following vector annihilation relation in the Fourier domain:
(5) |
Here for k = (k1, k2), 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 Λ1 and Λ2 is support limited to the dilation of Λ1 by Λ2, which we denote by Λ1 ∗ Λ2. Since is bandlimited to Λ1, when using samples of within Γ, we require Λ1 ∗ Λ2 = Γ; see Fig. 2.
The Fourier domain annihilation relations (5) can be compactly represented in matrix form as
(6) |
where , i = 1, 2, are matrices corresponding to the discrete 2-D convolution of , (omitting the irrelevant factor j2π) with a filter supported on Λ1, with output restricted to the index set Λ2. Here we use h to denote the vectorized version of a filter h[k], k ∈ Λ1. By our previous observations, the solutions h to (6) are given by the Fourier coefficients of a multiple of the minimal polynomial. Hence if the filter support Λ1 is larger than the minimal filter support Λ0, has a large nullspace and is low-rank. Specifically, in [8] we proved the following:
Proposition 1.
Suppose is built with filter size Λ1 ⊇ Λ0 satisfying Γ ⊇ 2Λ1 ∗ Λ0, then
where |Λ1| is the number of indices in Λ1 and |Λ1 : Λ0| is the number of integer shifts of Λ0 contained in Λ1.
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:
(7) |
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 k ∈ Γ, where
(8) |
for i = 1, 2. Here ωi(k) is the set of locations of the matrix containing copies of . Note that the set {Ak}k∈Γ 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:
(9) |
(10) |
(11) |
The constants in (11) are chosen so that . Using these definitions, we rewrite (7) as
(12) |
Several authors have shown that the performance of low-rank matrix recovery by nuclear norm minimization is dependent on the incoherence of the sampling basis with respect to the matrix to be to be recovered [2, 11]. Towards this end, we introduce two incoherence measures associated with defined solely in terms of the edge-set polynomial μ0. In the following, we set to be the 2-D Dirichlet kernel supported on Λ1, i.e., the function such that for all k ∈ Λ1 and zero otherwise. For any collection of N points P = {r1, …, rN} ⊂ [0, 1]2, we define the N × N Gram matrix G(P) by .
Definition 2.
Let μ0 be a trigonometric polynomial bandlimited to Λ0 (see (3))), and set R = |Λ1| − |Λ1 : Λ0|. Define the incoherence measure ρ1 by
(13) |
where σmin[G(P)] is the minimum singular value of G(P).
Put in words, among all possible arrangements of R points along the edge-set {μ0 = 0}, we seek the arrangement that gives the best conditioning of the matrix G(P), and call the resulting condition number ρ1. Intuitively, the optimal arrangement will maximize the minimum separation distance among the R points, and ρ1 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 ρ1.
Additionally, our results rely on another incoherence measure related to properties of the gradient of μ0:
Definition 3.
Let μ0 be a trigonometric polynomial bandlimited to Λ0 (see (3)). Normalize μ0 by . Define the incoherence measure ρ2 by
(14) |
where denotes the space of all trigonometric polynomials bandlimited to Λ1, and .
Note that ρ2 will be large when has several zeros, or equivalently, when μ0 has several critical points. Since μ0 must contain a critical point in every region defined by the complement of the edge-set, ρ2 will be large when the image has several distinct regions.
Now we state our main result:
Theorem 4.
Let f be specified by (1), whose edge-set is described by the zero-set of the trigonometric polynomial μ0 bandlimited to Λ0 (see (3)) with associated incoherence measures ρ1 and ρ2. Let Θ ⊂ Γ be an index set drawn uniformly at random within Γ. Then there exists a universal constant c > 0 such that the solution to (12) is exact with probability exceeding 1 − |Γ|−2, provided
(15) |
where R = rank and cs = |Γ|/|Λ1|.
Following the approach in [2], we can prove this result by constructing an approximate dual certificate using the well-known golfing scheme of [11]. The adaptation of the proof in [2] to the measurement operator (11) is straightforward, and these details are omitted for brevity. The essential difference of the above result and [2] is the characterization of the incoherency measures. The approach in [2] 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-R 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):
Proposition 5.
Consider of rank R corresponding to a piecewise constant function f whose edge set coincides with the zero set of μ0, let ρ1 and ρ2 be the incoherency measures of μ0, and set cs = |Γ|/|Λ1|. Then we have
(16) |
The proof relies on the following basis representations for the row and column spaces of :
Lemma 6.
Choose any N ≥ |Λ1| + |Λ0| − |Λ1 : Λ0| points {r1, …, rN} ⊂ {μ0 = 0}, and define B to be the collection of filters {h1, …, hN} where each are the Fourier coefficients of the translated Dirichlet kernel . Then there exists a subset of elements from B that is a basis for the row space of .
Lemma 7.
Let be any basis of the row space of , and set . Then is a basis of the column space of .
3. EXPERIMENTS
In Fig. 3 we demonstrate the recovery of a synthetic piecewise constant phantom [12] (181 × 181 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 [6]. The filter size was set to 33 × 25. 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 (255×255 sampling grid with 4 receiver coils, coil-compressed to a single channel) from 2-fold uniform random undersampling, using a filter size of 45 × 45. 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 [10] 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 [2] 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.
Acknowledgments
This work is supported by grants NIH 1R01EB019961-01A1, NSF CCF-1116067, ONR N00014-13-1-0202, and ACS RSG-11-267-01-CCE.
5. REFERENCES
- [1].Vetterli M, Marziliano P, and Blu T, “Sampling signals with finite rate of innovation,” Signal Processing, IEEE Transactions on, vol. 50, no. 6, pp. 1417–1428, 2002. [Google Scholar]
- [2].Chen Y and Chi Y, “Robust spectral compressed sensing via structured matrix completion,” Information Theory, IEEE Transactions on, vol. 60, no. 10, pp. 6576–6601, 2014. [Google Scholar]
- [3].Jin KH, Lee D, and Ye JC, “A general framework for compressed sensing and parallel MRI using annihilating filter based low-rank hankel matrix,” arXiv preprint arXiv:1504.00532, 2015. [Google Scholar]
- [4].Jin KH, Lee D, and Ye JC, “A novel k-space annihilating filter method for unification between compressed sensing and parallel mri,” in International Symposium on Biomedical Imaging (ISBI), April 2015. [Google Scholar]
- [5].Haldar JP, “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. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [6].Ongie G and Jacob M, “Recovery of piecewise smooth images from few fourier samples,” in Sampling Theory and Applications (SampTA), May 2015, pp. 543–547. [Google Scholar]
- [7].Ongie G and Jacob M, “Super-resolution MRI using finite rate of innovation curves,” IEEE International Symposium on Biomedical Imaging (ISBI), April 2015. [Google Scholar]
- [8].Ongie G and Jacob M, “Off-the-grid recovery of piecewise constant images from few fourier samples,” in arXiv:1510.00384, May 2015, pp. 543–547. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [9].Pan H, Blu T, and Dragotti PL, “Sampling curves with finite rate of innovation,” Signal Processing, IEEE Transactions on, vol. 62, no. 2, 2014. [Google Scholar]
- [10].Ongie G and Jacob M, “A fast algorithm for structured low-rank matrix recovery with applications to undersampled MRI reconstruction,” IEEE International Symposium on Biomedical Imaging (ISBI), April 2016. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [11].Gross D, “Recovering low-rank matrices from few coefficients in any basis,” Information Theory, IEEE Transactions on, vol. 57, no. 3, pp. 1548–1566, 2011. [Google Scholar]
- [12].Guerquin-Kern M, Lejeune L, Pruessmann KP, and Unser M, “Realistic analytical phantoms for parallel magnetic resonance imaging,” Medical Imaging, IEEE Transactions on, vol. 31, no. 3, pp. 626–636, 2012. [DOI] [PubMed] [Google Scholar]