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.

## 2 Theory

### 2.1 Signal Model: 2-D Piecewise Constant Images

We consider the recovery of a piecewise constant function

 f(r)=N∑i=1ai χΩi(r),   for all   r=(x,y)∈[0,1]2, (1)

where , and denotes the characteristic function of the set . We assume the Fourier samples specified by

 ^f[k]=∫[0,1]2f(r)e−j2πk⋅r; k∈Z2, (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:

 μ0(r)=∑k∈Λ0c[k]ej2πk⋅r,∀r∈[0,1]2, (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

 μ∇f=0, (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:

 ∑k∈Λ1ˆ∇f[ℓ−k]ˆμ[k]=0,  ∀ ℓ∈Λ2. (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

 T(^f)=[T1(^f)T2(^f)]h=0, (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

 R:=rank T(ˆf)=|Λ1|−|Λ1|Λ0|

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:

 minimizeˆg ∥T(ˆg)∥∗ % subject to ˆg[k]=ˆf[k];k∈Θ (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

 (Ai,k)α,β = ⎧⎨⎩ki|k|√|ωi(k)|, if (α,β)∈ωi(k)0else (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:

 AΘ(X) =∑k∈Θ⟨Ak,X⟩Ak (9) A⊥(X) =I−∑k∈Γ⟨Ak,X⟩Ak (10) QΘ =|Γ||Θ|AΘ+A⊥ (11)

The constants in (11) are chosen so that . Using these definitions, we rewrite (7) as

 minimizeX ∥X∥∗ subject to QΘ(X)=QΘ(T(^f)) (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

 1ρ1=maxP⊂{μ0=0}|P|=Rσmin[G(P)] (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

 1ρ2=minγ∈BΛ1∥γ∥2=1∫[0,1]2|γ(r)|2|∇μ0(r)|2dr∥ˆ∇μ0∥21 (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

 |Θ|>cρ1ρ2Rcslog4|Γ|. (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 T(^f) 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

 maxk∈Γ{∥PUAk∥,∥PVAk}∥2F≤ρ1ρ2Rcs|Γ| (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 .