Structured low-rank recovery of piecewise constant signals with performance guarantees

Structured low-rank recovery of piecewise constant signals with performance guarantees

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.

Structured low-rank recovery of piecewise constant signals with performance guarantees

Greg Ongie, Sampurna Biswas, and Mathews Jacobthanks: 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.

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, 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 [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.

Fig. 1: Annihilation of a piecewise constant function in the spatial (top) and Fourier (bottom) domain.

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 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:

(3)

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 [8] that when is any bandlimited function that vanishes on the edgeset , the gradient satisfies the property

(4)

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 [8]. The spatial domain relation (4) translates directly to the following vector annihilation relation in the Fourier domain:

(5)

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

(6)

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 [8] we proved the following:

Proposition 1.

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 .

Fig. 2: Fourier domain support sets used in constructing the structured matrix . The grid represent a rectangular sampling window in (index marked in black). (left) is a collection of uniform random sampling locations. (right) is the Fourier support of the minimal annihilating polynomial, is the assumed filter size used in constructing , and is the set of valid convolutions that satisfies . Note that has dimensions .

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 , where

(8)

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:

(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 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 .

Definition 2.

Let be a trigonometric polynomial bandlimited to (see (3)), and set . Define the incoherence measure by

(13)

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 :.

Definition 3.

Let be a trigonometric polynomial bandlimited to (see (3)). Normalize such that . Define the incoherence measure by

(14)

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:

Theorem 4.

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

(15)

where and .

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- 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 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

(16)

The proof relies on the following basis representations for the row and column spaces of :

Lemma 6.

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 .

Lemma 7.

Let be any basis of the row space of , and set . Then is a basis of the column space of .

Fully sampled

TV regularized recovery

Structured low-rank recovery
Fig. 3: Recovery of synthetic MRI data from 20-fold variable density undersampling.

Fully sampled (top) Subsampled -space (bottom)

TV regularized recovery SNR=17.8dB

Structured low-rank recovery SNR=19.0dB

Fig. 4: Recovery of real MRI data from 2-fold random uniform undersampling. Error images shown below.

3 Experiments

In Fig. 3 we demonstrate the recovery of a synthetic piecewise constant phantom [12] ( 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 . 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 [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.

References

  • [1] 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.
  • [2] 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.
  • [3] 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.
  • [4] ——, “A novel k-space annihilating filter method for unification between compressed sensing and parallel mri,” in International Symposium on Biomedical Imaging (ISBI), April 2015.
  • [5] 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.
  • [6] 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.
  • [7] ——, “Super-resolution MRI using finite rate of innovation curves,” IEEE International Symposium on Biomedical Imaging (ISBI), April 2015.
  • [8] ——, “Off-the-grid recovery of piecewise constant images from few fourier samples,” in arXiv:1510.00384, May 2015, pp. 543–547.
  • [9] 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.
  • [10] 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.
  • [11] 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.
  • [12] 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.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
49341
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description