Convex recovery of continuous domain piecewise constant images from non-uniform Fourier samples

TL;DR

This paper introduces a convex matrix completion method for recovering continuous domain piecewise constant images from non-uniform Fourier samples, leveraging low-rank structures induced by edge set properties.

Contribution

It generalizes super-resolution techniques to piecewise constant images, modeling edges as zero levelsets of bandlimited functions and formulating a structured low-rank matrix completion problem.

Findings

Exact recovery is probable under an incoherency condition.

The incoherency depends on the edge set geometry, affecting sampling requirements.

The method extends super-resolution to more complex image structures.

Abstract

We consider the recovery of a continuous domain piecewise constant image from its non-uniform Fourier samples using a convex matrix completion algorithm. We assume the discontinuities/edges of the image are localized to the zero levelset of a bandlimited function. This assumption induces linear dependencies between the Fourier coefficients of the image, which results in a two-fold block Toeplitz matrix constructed from the Fourier coefficients being low-rank. The proposed algorithm reformulates the recovery of the unknown Fourier coefficients as a structured low-rank matrix completion problem, where the nuclear norm of the matrix is minimized subject to structure and data constraints. We show that exact recovery is possible with high probability when the edge set of the image satisfies an incoherency property. We also show that the incoherency property is dependent on the geometry of the edge set curve, implying higher sampling burden for smaller curves. This paper generalizes recent work on the super-resolution recovery of isolated Diracs or signals with finite rate of innovation to the recovery of piecewise constant images.

Tables

Table 1. TABLE I: Comparison of proposed scheme with discrete total variation minimization

	TV-minimization	Proposed
Spatial domain	discrete	continuous
Derivative operator	finite differences	exact derivative
Singularity set	discrete points	connected curves
Frequency domain	discrete	discrete
Frequency weighting $w_i\left[𝒌\right]$	$1-e^{j2\pi k_i/N_i}$	$j2\pi k_i$
Lifted matrix structure	two-level circulant	two-level Toeplitz
Rank of lifted matrix	sparsity of	bandwidth of
	discrete gradient	edge set

Equations

$$f(\mathbf{r})=\sum _{i=1}^N{a}_i{1}_{U_i}(\mathbf{r}),\text{ for all }\mathbf{r}=(x,y)\in {[0,1]}^2,$$

$$\hat{f}[\mathbf{k}]={\int }_{[0,1{]}^2}f(\mathbf{r})e^{-j2\pi \mathbf{k}\cdot \mathbf{r}};\mathbf{k}\in \Omega \subset {\mathbb{Z}}^2.$$

$$E=\{\mathbf{r}\in [0,1{]}^2:\mu (\mathbf{r})=0\},\text{with}\mu (\mathbf{r})=\sum _{\mathbf{k}\in \Lambda }c[\mathbf{k}]e^{j2\pi \mathbf{k}\cdot \mathbf{r}},$$

$$\mu \nabla f=0,$$

$$\sum _{\mathbf{k}\in \Lambda }\widehat{\nabla f}[\mathbf{\ell }-\mathbf{k}]\hat{\mu }[\mathbf{k}]=0,\forall \mathbf{\ell }\in {\mathbb{Z}}^2.$$

$$\mathcal{T}(\hat{f})\mathbf{h}=\left[\begin{array}{c}{\mathcal{T}}_x(\hat{f})\\ {\mathcal{T}}_y(\hat{f})\end{array}\right]\mathbf{h}=0,$$

$$\mathbf{h}=(({\hat{\mu }}_0\ast \hat{\gamma })[\mathbf{k}]:\mathbf{k}\in {\Lambda }_1)$$

$$\text{rank}\mathcal{T}(\hat{f})\le |{\Lambda }_1|-|{\Lambda }_1:{\Lambda }_0|$$

$${\mu }_0\nabla g=0\text{subject to}\hat{g}[\mathbf{k}]=\hat{f}[\mathbf{k}]\text{for all}\mathbf{k}\in \Gamma ,$$

$$\underset{\hat{g}[\mathbf{k}],\mathbf{k}\in \Gamma }{\min }\Vert \mathcal{T}(\hat{g}){\Vert }_{\ast }\text{mbox}subjectto\hat{g}[\mathbf{k}]=\hat{f}[\mathbf{k}]\text{for all}\mathbf{k}\in \Omega$$

$$\mathcal{T}(\hat{f})=\sum _{\mathbf{k}\in \Gamma /\{0\}}\hat{f}[\mathbf{k}]w[\mathbf{k}]{\mathbf{A}}_{\mathbf{k}}$$

$$\underset{\mathbf{k}\in \Gamma }{\max }\Vert {\mathcal{P}}_U{\mathbf{A}}_{\mathbf{k}}{\Vert }_F^2$$

$$\underset{\mathbf{k}\in \Gamma }{\max }\Vert {\mathcal{P}}_V{\mathbf{A}}_{\mathbf{k}}{\Vert }_F^2$$

$$|\Omega |>c\rho c_{s}R{\log }^4|\Gamma |,$$

$$\left(\frac{|{\Lambda }_1|-|{\Lambda }_1:{\Lambda }_0|}{|{\Lambda }_1|}\right)|\Gamma |\le \left(\frac{{\alpha }^2-(\alpha -1{)}^2}{{\alpha }^2}\right)|\Gamma |\le \frac{2|\Gamma |}{\alpha }.$$

$$c_{s}R=O\left(\alpha |{\Lambda }_0|\right).$$

$${\hat{f}}_n[\mathbf{k}]=\hat{f}[\mathbf{k}]+\eta [\mathbf{k}],\mathbf{k}\in \Omega ,$$

$$\underset{\hat{g}}{\min }\Vert \mathcal{T}(\hat{g}){\Vert }_{\ast }\text{subject to}\Vert {\mathcal{P}}_{\Omega }(\widehat{f_n}-\hat{g}){\Vert }_2\le \delta .$$

$$\Vert \mathcal{T}(\hat{f})-\mathcal{T}(\hat{g}){\Vert }_F\le 5|\Gamma {|}^2\delta .$$

$$\Vert f_n-f{\Vert }_{L^2}^2={\left({\int }_{[0,1{]}^2}|f_n(\mathbf{r})-f(\mathbf{r}){|}^{2}d\mathbf{r}\right)}^{\frac{1}{2}}\le \delta .$$

$${\mathbf{E}}_{\text{row}}(P):=\frac{1}{\sqrt{|{\Lambda }_1}|}\left(\begin{array}{ccc}{e}^{j2\pi {\mathbf{k}}_1\cdot {\mathbf{r}}_1}& \dots & e^{j2\pi {\mathbf{k}}_1\cdot {\mathbf{r}}_R}\\ ⋮& & ⋮\\ e^{j2\pi {\mathbf{k}}_{|{\Lambda }_1|}\cdot {\mathbf{r}}_1}& \dots & e^{j2\pi {\mathbf{k}}_{|{\Lambda }_1|}\cdot {\mathbf{r}}_R}\end{array}\right)$$

$${\mathbf{E}}_{\text{col}}(P)=\frac{1}{\sqrt{|{\Lambda }_2}|}\left(\begin{array}{ccc}\frac{w_{1,x}}{\Vert {\mathbf{w}}_1\Vert }{e}^{j2\pi {\mathbf{k}}_1\cdot {\mathbf{r}}_1}& \dots & \frac{w_{R,x}}{\Vert {\mathbf{w}}_R\Vert }{e}^{j2\pi {\mathbf{k}}_1\cdot {\mathbf{r}}_R}\\ ⋮& & ⋮\\ \frac{w_{1,x}}{\Vert {\mathbf{w}}_1\Vert }{e}^{j2\pi {\mathbf{k}}_{|{\Lambda }_2|}\cdot {\mathbf{r}}_1}& \dots & \frac{w_{R,x}}{\Vert {\mathbf{w}}_R\Vert }{e}^{j2\pi {\mathbf{k}}_{|{\Lambda }_2|}\cdot {\mathbf{r}}_R}\\ & & \\ \frac{w_{1,y}}{\Vert {\mathbf{w}}_1\Vert }{e}^{j2\pi {\mathbf{k}}_1\cdot {\mathbf{r}}_1}& \dots & \frac{w_{R,y}}{\Vert {\mathbf{w}}_R\Vert }{e}^{j2\pi {\mathbf{k}}_1\cdot {\mathbf{r}}_R}\\ ⋮& & ⋮\\ \frac{w_{1,y}}{\Vert {\mathbf{w}}_1\Vert }{e}^{j2\pi {\mathbf{k}}_{|{\Lambda }_2|}\cdot {\mathbf{r}}_1}& \dots & \frac{w_{R,y}}{\Vert {\mathbf{w}}_R\Vert }{e}^{j2\pi {\mathbf{k}}_{|{\Lambda }_2|}\cdot {\mathbf{r}}_R}\end{array}\right),$$

$$(\mathbf{G}(P){)}_{i,j}=\frac{1}{|{\Lambda }_1|}{D}_{{\Lambda }_1}({\mathbf{r}}_i-{\mathbf{r}}_j),1\le i,j\le R,$$

$$\rho =\underset{\begin{array}{c}P\subset \{{\mu }_0=0\}\\ |P|=R\end{array}}{\min }\frac{1}{{\lambda }_{min}[\mathbf{G}(P)]},$$

$$\rho \le {\left(1-\frac{1}{\sqrt{|{\Lambda }_1|}\Delta }\right)}^{-2},$$

$$\underset{\mathbf{u}\in {\mathbb{C}}^N}{\min }TV(\mathbf{u})\text{subject to}{P}_{\Omega }(\mathbf{F}\mathbf{u})=P_{\Omega }(\mathbf{F}{\mathbf{u}}_0).$$

$$TV(\mathbf{u})=\sum _{i=1}^N(|({\partial }_1\mathbf{u}{)}_i{|}^2+|({\partial }_2\mathbf{u}{)}_i{|}^2{)}^{\frac{1}{2}}$$

$$\underset{\mathbf{v}}{\min }\Vert \mathcal{C}(\mathbf{F}\mathbf{u}){\Vert }_{\ast }\text{subject to}{P}_{\Omega }(\mathbf{F}\mathbf{u})=P_{\Omega }(\mathbf{F}{\mathbf{u}}_0).$$

$$\mathcal{C}(\mathbf{F}\mathbf{u})=\left[\begin{array}{c}{\mathcal{C}}_x(\mathbf{F}\mathbf{u})\\ {\mathcal{C}}_y(\mathbf{F}\mathbf{u})\end{array}\right]\in {\mathbb{C}}^{2N\times N}$$

$$\mathcal{S}(f)\phi =D_{{\Lambda }_2}\ast (\phi \nabla f)\in B_{{\Lambda }_1}^2,$$

$$\mathcal{S}(f{)}^{\ast }\mathbf{\rho }=D_{{\Lambda }_1}\ast (\mathbf{\rho }\cdot \nabla f)\in B_{{\Lambda }_1}$$

Full text

I Introduction

The direct recovery of continuous domain signals by convex optimization is emerging as a powerful alternative to traditional discrete domain compressed sensing [1, 2, 3]. The ability of these continuous domain “off-the-grid” schemes to minimize discretization errors makes them attractive in practical applications, where only the low-pass measurements of the signal are available. The history of such continuous domain signal recovery algorithms dates back to Prony [4], where the recovery of a linear combination of exponentials from uniform samples is considered. Prony-like algorithms recover the signal by estimating an annihilating polynomial whose zeros correspond to the frequencies of the exponentials. Work by Liang et al. [5, 6] and the finite rate of innovation (FRI) framework [7] extended Prony-like 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 [3, 8, 9, 10, 11] by exploiting the low-rank structure of an enhanced matrix built from Fourier data (e.g., a Hankel matrix in 1-D). Recovery guarantees exists for certain classes of these signals when the singularities are isolated and well-separated [2, 3, 12].

The signal models discussed above have limited flexibility in exploiting the extensive additional structure present in multidimensional imaging problems. In particular, the edges in multidimensional images are connected and can be modeled as smooth curves or surfaces. While discrete image representations to capture this structure have been the subject extensive research [13, 14], similar continuous domain representations have attracted less attention. We recently introduced a novel framework recover piecewise polynomial images, whose edges are localized to smooth curves, from their uniform [15, 16] and non-uniform [11] Fourier samples; our framework generalizes a recent extension of FRI models to curves [17]. We assume that the partial derivatives of the signal vanish outside the zero level-set of a bandlimited function, which is only true for piecewise smooth signals. This relation translates to a linear system of convolution equations 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. We have introduced theoretical guarantees for the recovery of such images from uniform samples [16, 15]. Our earlier work has shown that the structured matrix built from the Fourier coefficients of piecewise constant images is low-rank [11, 16], which we used to recover the image from its non-uniform Fourier samples with good performance in practical applications. We have also introduced an computationally efficient algorithm termed as GIRAF, which works on the original signal samples rather than the structured high-dimensional matrix [18, 19]; the computational complexity of this algorithm is comparable to discrete total variation regularization, which makes this scheme readily applicable to large-scale imaging problems, such as undersampled dynamic magnetic resonance image reconstruction [20].

The main focus of the present paper is to introduce theoretical guarantees on the recovery of continuous domain piecewise constant images from non-uniform Fourier samples via a convex structured low-rank matrix completion algorithm. Our main result shows number of non-uniform samples to recovery the image is proportional to the complexity of the edge set, as measured by the bandwidth of the edge set function, and an incoherence measure related to the edge set geometry. We additionally show that the recovery is robust to noise and model-mismatch.

The proof of the main result builds off of [3], which proved similar recovery guarantees for the recovery of multi-dimensional isolated Diracs from non-uniform Fourier samples by minimizing the nuclear norm of an “enhanced” multi-level Hankel matrix. This work showed that the number of samples necessary for recovery depends the number of Diracs and on an incoherence measure of the signal, that can be defined solely in terms of the relative locations of the Diracs. However, the theory in [3] relies heavily on an 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 singularities are isolated and finite. Since the singularities in the proposed class of piecewise constant images (i.e., the image edges) are not isolated nor finite, the recovery guarantees in [3] cannot be directly extended to our setting. Instead, to achieve our result, we give a new characterization of the row and column spaces of the structured matrix arising in our setting. We show this new characterization allows us to derive an incoherence measure based solely on geometric properties of the edge set. In particular, we derive an upper bound for the incoherence measure that is related to the size of edge set curve. The results show that high sampling burden is associated with the estimation of images with smaller piecewise constant regions, which is consistent with intuition.

We note that the signal models in [3, 1, 2] do not include the class of piecewise constant images considered in this work. In particular, all of the above models assume the discontinuities to be finite in number and well separated, unlike in our setting. Recently, [12] adapted the results in [3] to introduce recovery guarantees for Fourier interpolation of a variety of finite-rate-of-innovation signal models [7], including piecewise constant functions. However, these results are limited to the 1-D setting and share the assumption than the discontinuities/innovations of the signal are finite and isolated. Furthermore, the structured matrix lifting considered in this work is different than those considered in [3] and [12]. Specifically, the structured matrix lifting in this work consists of two vertically concatenated multi-level Toeplitz matrices (i.e., block Toeplitz with Toeplitz blocks), whose entries are built from the weighted Fourier coefficients of the images. This is substantially different from the structured matrix liftings considered in [3] (unweighted, one block, single block multilevel Hankel) and [12] (weighted, one block, single-level Hankel). Finally, we note that a preliminary version of the results presented in this has been published previously in the conference paper [21] without proofs. The present work includes considerably more details and proofs, and major improvements to the main theorem.

I-A Notation

Bold lower-case letters $\bm{x}$ are used to indicate vector quantities, bold upper-case $\bm{X}$ to denote matrices, and calligraphic script $\mathcal{X}$ for general linear operators. We typically reserve lower-case greek letters $\mu,\gamma$, etc. for trigonometric polynomials (3) and upper-case greek letters $\Lambda,\Omega,$ etc. for their coefficient index sets, i.e. finite subsets of the integer lattice $\mathbb{Z}^{2}$, with cardinality denoted by $|\Lambda|$. We write $\Lambda+\Omega$ for the dilation of the index set $\Omega$ by $\Lambda$, i.e. the Minkowski sum $\{\bm{k}+\bm{\ell}:\bm{k}\in\Lambda,~{}\bm{\ell}\in\Omega\}$, and write $2\Lambda$ to mean $\Lambda+\Lambda$, $3\Lambda=2\Lambda+\Lambda$, etc. We also denote the contraction of $\Omega$ by $\Lambda$ by $\Omega\,{:}\,\Lambda=\{\bm{\ell}\in\Omega\,{:}\,\bm{\ell}-\bm{k}\in\Omega\text{ for all }\bm{k}\in\Lambda\}$.

II Background

II-A 2-D Piecewise Constant Images with Bandlimited Edges

In this work we consider a continuous domain piecewise constant model for images,

[TABLE]

where $a_{i}\in\mathbb{C}$, $1_{U}$ denotes the characteristic function of the set $U$, and each $U_{i}\subset[0,1]^{2}$ is a simply connected regions with piecewise smooth boundaries $\partial U_{i}$. We study the recovery of such an image from a sampling of its Fourier coefficients $\hat{f}$ specified by

[TABLE]

Following [16], we further assume that the edge set of the piecewise constant image, specified by $E:=\cup_{i}\partial U_{i}$, coincides with the zero set of a 2-D bandlimited function:

[TABLE]

where the coefficients $c[\bm{k}]\in\mathbb{C}$, and ${\Lambda}$ is a finite subset of $\mathbb{Z}^{2}$. We call any function $\mu$ in the form (3) a trigonometric polynomial, and we say $\mu$ is bandlimited to $\Lambda$, i.e., the Fourier coefficients $\hat{\mu}$ are supported within $\Lambda$. For short, we will write $\{\mu=0\}$ for the zero set of $\mu$ considered as a subset of $[0,1]^{2}$.

Define the degree of a trigonometric polynomial $\mu$, denoted by $deg(\mu)=(K,L)$ to be the linear dimensions of the smallest rectangle containing the support set $\{\bm{k}:\hat{\mu}[\bm{k}]\neq 0\}$. In [16] we proved that for every curve $E$ given by the zero set of a trigonometric polynomial, there exists a unique minimal degree trigonometric polynomial111More precisely, $\mu_{0}$ is unique up to multiplication by a phase factor $e^{j2\pi\bm{k}\cdot\bm{r}}$ for some $\bm{k}\in\mathbb{Z}^{2}$. $\mu_{0}$ such that $E=\{\mu_{0}=0\}$ and if $\mu$ is any other trigonometric polynomial with $\{\mu_{0}=0\}\subset\{\mu=0\}$, then $deg(\mu_{0})\leq deg(\mu)$ entrywise. By extension, we define the degree of a curve $E$ to be equal to the degree of of its minimial degree polynomial $\mu_{0}$. We also say the curve $E$ is bandlimited to $\Lambda_{0}\subset\mathbb{Z}^{2}$, where $\Lambda_{0}$ is the minimal rectangular index set containing the support of $\hat{\mu}$. Intuitively, the degree/bandwidth of a curve gives a quantitative measure of its complexity. For example, in [16] we show the number of connected components of a curve is bounded in terms of its degree.

II-B Recovery from uniform Fourier samples

We have shown in [16] that when $\mu$ is any trigonometric polynomial that vanishes on the edge set of the piecewise constant image $f$, the gradient $\bm{\nabla}f=\left(\partial_{x}f,\partial_{y}f\right)$ satisfies the property

[TABLE]

where equality in (4) is understood in the sense of distributions (see, e.g., [22]). The spatial domain annihilation relation (4) translates directly to the following convolution annihilation relation in Fourier domain:

[TABLE]

Here $\widehat{\bm{\nabla}f}[\bm{k}]=j2\pi(k_{x}\widehat{f}[\bm{k}],k_{y}\widehat{f}[\bm{k}])$ for $\bm{k}=(k_{x},k_{y})$. Note the equations in (5) are linear with respect to the coefficients $\widehat{\mu}$.

Suppose we have access to samples of the Fourier coefficients $\hat{f}$ on a finite rectangular grid $\Gamma\subset\mathbb{Z}^{2}$, and suppose $\mu$ is bandlimited to $\Lambda_{1}\subset\Gamma$. Then we can build the system of equations in (5) for all $\bm{\ell}$ belonging to the index set $\Lambda_{2}\subset\Gamma$, where $\Lambda_{2}$ is the set of all integer shifts of $\Lambda_{1}$ contained in $\Gamma$. In this case (5) can be compactly represented in matrix form as

[TABLE]

where $\mathcal{T}_{x}(\hat{f}\,),\mathcal{T}_{y}(\hat{f}\,)\in\mathbb{C}^{|\Lambda_{2}|\times|\Lambda_{1}|}$ are matrices corresponding to the discrete 2-D convolution with the arrays $k_{x}\widehat{f}[k_{x},k_{y}]$ and $k_{x}\widehat{f}[k_{x},k_{y}]$ for $(k_{x},k_{y})\in\Gamma$, respectively (after omitting the inconsequential factor $j2\pi$). Here we use $\bm{h}$ to denote the vectorized version of the filter $(\hat{\mu}[\bm{k}]:\bm{k}\in\Lambda_{1})$, where the index set $\Lambda_{1}$ is called the filter support. The matrices $\mathcal{T}_{x}(\hat{f}\,)$ and $\mathcal{T}_{y}(\hat{f}\,)$ have a block Toeplitz with Toeplitz blocks structure. See Figure 2 for an illustration of the construction of $\mathcal{T}(\hat{f}\,)$.

Equation (6) shows that $\mathcal{T}(\hat{f}\,)$ is rank deficient, since it has the non-trivial vector $\bm{h}$ in its nullspace. In addition, when the filter support $\Lambda_{1}$ defining $\mathcal{T}(\hat{f}\,)$ is sufficiently big, we can also show $\mathcal{T}(\hat{f}\,)$ is low-rank. This is because if $\mu_{0}$ is the minimal degree polynomial for the edge set, then any multiple of $\mu=\gamma\cdot\mu_{0}$ bandlimited to $\Lambda_{1}$ will satisfy the annihilation equation (4). In Fourier domain, this means the vector

[TABLE]

is in the nullspace of $\mathcal{T}(\hat{f}\,)$. Hence if the filter support $\Lambda_{1}$ is larger than support $\Lambda_{0}$ of $\mu_{0}$, $\mathcal{T}(\hat{f}\,)$ has a large nullspace and is low-rank. The following result from [16] gives an exact characterization of the rank of $\mathcal{T}(\hat{f}\,)$, which will be important for this work:

Theorem 1**.**

[16]** Suppose $f$ is a piecewise constant image (1) whose edge set $E=\{\mu_{0}=0\}$ is the zero set of a trigonometric polynomial $\mu_{0}$ bandlimited to $\Lambda_{0}$. Let $\mathcal{T}(\hat{f}\,)$ be built with filter size $\Lambda_{1}\supseteq\Lambda_{0}$, then

[TABLE]

where $|\Lambda_{1}|$ is the number of indices in $\Lambda_{1}$ and $|\Lambda_{1}\colon\Lambda_{0}|$ is the number of integer shifts of $\Lambda_{0}$ contained in $\Lambda_{1}$. Moreover, equality holds in (8) if $\Gamma\supseteq 2\Lambda_{1}+\Lambda_{0}$ and if the edge set does not contain any singular points. In this case, the nullspace of $\mathcal{T}(\hat{f}\,)$ consists of all vectors in the form (7).

Note that $R:=|\Lambda_{1}|-|\Lambda_{1}\colon\Lambda_{0}|$ is a measure of the bandwidth of $\mu_{0}$ and hence is indicative of the complexity of the edge set curve $E=\{\mu_{0}=0\}$. In the remainder of this work we assume the conditions in Theorem 1 that guarantee the equality $\text{rank}~{}\mathcal{T}(\hat{f}\,)=R$ holds, in particular $\Gamma\supseteq 2\Lambda_{1}+\Lambda_{0}$.

If we take $\Lambda_{1}=\Lambda_{0}$, the above result shows Fourier samples of $\hat{f}$ in $\Gamma\supseteq 3\Lambda_{0}$ is sufficient for the recovery of the minimal degree polynomial $\mu_{0}$, since in this case $\hat{\mu}_{0}$ can be identified as the unique non-trivial nullspace vector of $\mathcal{T}(\hat{f}\,)$. The following theorem states that once $\mu_{0}$ is available, $f$ is the unique solution to the annihilation equations (4) and (5):

Theorem 2**.**

[16*]**.

Suppose $f$ is a piecewise constant image (1) whose edge set $E=\{\mu_{0}=0\}$ is the zero set of a trigonometric polynomial $\mu_{0}$ bandlimited to $\Lambda_{0}$. Suppose the Fourier sampling set $\Gamma\supseteq\Lambda_{0}$. If $g\in L^{1}([0,1]^{2})$ satisfies*

[TABLE]

then $g=f$ almost everywhere.

In principle, this result allows us to solve for the amplitudes of regions of the piecewise constant function $f$ by plugging in the known $\mu_{0}$ into the equation (9) and solving a linear system, similar to Prony’s method. However, for complicated piecewise constant images with many regions, it may be more practical to use the approximations introduced in [16].

III Recovery from non-uniform Fourier samples

The theory presented in Section I shows that the exact recovery of a continuous domain piecewise constant image with a bandlimited edge set is possible when we collect Fourier samples of the image on a sufficiently large uniform grid in Fourier domain. However, the recovery procedure breaks down when we have non-uniform or missing samples, which is often the case in practical settings, e.g., compressed sensing MRI [23]. Therefore, we propose and analyze a method to interpolate the missing samples to a uniform grid in Fourier domain, which guarantees full recovery of the image in spatial domain.

Recall that Theorem 1 says that the structured matrix $\mathcal{T}(\hat{f}\,)$ built from the Fourier coefficients $\hat{f}[\bm{k}],\bm{k}\in\Gamma$, where $\Gamma\subset\mathbb{Z}^{2}$ is a uniform rectangular grid, is known to be low-rank precisely when $f$ is a piecewise constant image with a bandlimited edge set. Hence we propose to recover $\hat{f}[\bm{k}],\bm{k}\in\Gamma$ from its samples at non-uniform locations $\Omega\subset\Gamma$ as the solution to the convex matrix completion problem:

[TABLE]

where $\|\cdot\|_{*}$ denotes the nuclear norm, i.e., the sum of the singular values of a matrix, which is the convex relation of the rank functional. Note that (10) is different than the standard low-rank matrix completion setting studied in [24, 25] in that the low-rank matrix $\mathcal{T}(\hat{f}\,)$ is structured and parameterized by the coefficient vector $\hat{f}$. Similar structured low-rank matrix completion schemes have been proposed for the recovery of signals from non-uniform Fourier samples [3, 12] and used with empirical success in MRI applications [26, 10, 11]. The main focus of this paper is to determine the sufficient number of samples that will ensure exact recovery of the Fourier coefficients of $f$ on the reconstruction grid $\Gamma$ with high probability.

III-A Role of incoherence

Several authors have shown that the sufficient number of samples for low-rank matrix recovery by nuclear norm minimization to succeed is dependent on the incoherence of the sampling basis with respect to the matrix to be to be recovered [25, 3]. Similarly, our results depend on an incoherence measure derived from the structure of the matrix $\mathcal{T}(\hat{f}\,)$ and properties of the piecewise constant image $f$. In particular, define $\mathcal{P}_{U}$ and $\mathcal{P}_{V}$ to be the orthogonal projections onto the column space and row space of $\mathcal{T}(\hat{f}\,)$, respectively, i.e., if $\mathcal{T}(\hat{f}\,)=\bm{U}\bm{\Sigma}\bm{V}^{*}$ is the rank-$R$ singular value decomposition then $\mathcal{P}_{U}\bm{X}=\bm{U}\bm{U}^{*}\bm{X}$, $\mathcal{P}_{V}\bm{X}=\bm{X}\bm{V}\bm{V}^{*}$. In Appendix B, we show that the structured matrix $\mathcal{T}(\hat{f}\,)$ can be expanded using orthonormal basis of matrices $\mathbf{A}_{k}$ such that

[TABLE]

where $w[\bm{k}],\bm{k}\in\Gamma/\{0\}$ are a set of positive weights that do not depend on $\hat{f}$. Similar to results in [25, 3, 12], we prove that nuclear norm minimization (10) recovers the exact low-rank solution with high probability provided we can uniformly bound the norms of the projections of the sampling basis matrices $\bm{A}_{\bm{k}}$ onto the row and column spaces of $\mathcal{T}(\hat{f}\,)$:

Proposition 3**.**

Consider $\mathcal{T}(\hat{f}\,)$ of rank $R$ corresponding to a piecewise constant function $f$ whose edge set coincides with the zero set of $\mu_{0}$, let $\rho$ be the incoherency measure associated with $\mu_{0}$ to be defined in the sequel, and set $c_{s}=|\Gamma|/|\Lambda_{1}|$. Then we have

[TABLE]

The proof in Section VIII-F relies on the row and column spaces of $\mathcal{T}(\hat{f}\,)$ derived in Lemma 8 and Lemma 6 in the next section. These results will be used in the derivation of the main theorem in Section IX.

III-B Main Results

We now present our main results, which determine the sufficient number of random Fourier samples for the convex structured low-rank matrix completion program (10) to succeed with high probability. Our first theorem addresses the case of recovery from noiseless Fourier samples:

Theorem 4**.**

Let $f$ be a continuous domain piecewise constant image (1), whose edge-set is described by the zero-set of the trigonometric polynomial $\mu_{0}$ bandlimited to $\Lambda_{0}$ (see (3)). Let $\Omega\subset\Gamma$ be an index set drawn uniformly at random within $\Gamma$. Then there exists a universal constant $c>0$ such that the solution to (10) is $\widehat{f}$ with probability exceeding $1-|\Gamma|^{-2}$, provided

[TABLE]

where $R=|\Lambda_{1}|-|\Lambda_{1}:\Lambda_{0}|=\text{rank}~{}\mathcal{T}(\hat{f}\,)$, $c_{s}=|\Gamma|/|\Lambda_{1}|$, $c$ is a universal constant, and $\rho\geq 1$ is an incoherence measure depending on the geometry of the edge-set, to be defined in the sequel.

To better understand the dependence of the bound in (14) on the filter size $\Lambda_{1}$ and the edge set bandwidth $\Lambda_{0}$, assume for simplicity that $\Lambda_{1}$ is some dilation of $\Lambda_{0}$, that is, $\Lambda_{1}=\alpha\Lambda_{0}$, where $\alpha>1$ is an integer. In this case, the factor $c_{s}\,R$ in (14) simplifies to

[TABLE]

Therefore, assuming the other constants in (14) are fixed, the number of measurements sufficient for exact recovery is proportional to the reciprocal of the dilation factor $\alpha$. This suggests taking the filter size $\Lambda_{1}$ to be as large as allowed by Theorem 4. Namely, $\Lambda_{1}$ should satisfy $2\Lambda_{1}+\Lambda_{0}=\Gamma$, i.e., the side-lengths of filter support $\Lambda_{1}$ should be roughly half those of the reconstruction grid $\Gamma$. Fixing the filter support $\Lambda_{1}$ to obey this bound, then $\Gamma=(2\alpha+1)\Lambda_{0}$, and so $|\Gamma|\leq(2\alpha+1)^{2}|\Lambda_{0}|$. Inserting this bound into (15) gives

[TABLE]

Combined with (14), this shows that the number of measurements sufficient for exact recovery is on the order of $|\Lambda_{0}|$, up to incoherence and log factors.

The proof of Theorem 4, detailed in Appendix B, is in line with the approach of [3]. In particular, we prove the result by constructing an approximate dual certificate using the well-known “golfing scheme” of [25]. The main differences between in the proof of the above result and that in [3] results from the differences in the matrix structure and hence the characterization of the incoherency between the row and column subspaces of $\mathcal{T}(\hat{f}\,)$ with the sampling basis. In particular, the matrix $\mathcal{T}(\hat{f}\,)$ we consider is obtained by stacking two block Toeplitz with Toeplitz blocks (BTTB) matrices whose entries are the weighted Fourier coefficients of $f$, as opposed to a single unweighted BTTB matrix in [3]. The approach in [3] relies on an explicit low-rank factorization of a BTTB matrix in terms of Vandermonde-like matrices222The structured matrices considered in [3] are block Hankel with Hankel block matrices (BHHB), but this difference is purely cosmetic: every BTTB matrix can be re-expressed as BHHB after a permutation of its rows and columns. In particular, the Vandermonde-like factorization of BHHB matrices in [3] carries over to BTTB matrices.. Since this factorization is not available in our setting, we use algebraic properties of trigonometric polynomials to give a new characterization of the row and column spaces of the matrix. In particular, we show in Section IV that similar Vandermonde-like basis matrices exist for the row and column space of the lifted matrix, and use these to derive a related incoherence measure that satisfies the bounds in Prop. 3.

III-C Recovery in the presence of noise and model-mismatch

We now generalize (68) to the setting where we have noisy or corrupted Fourier samples

[TABLE]

where $\eta[\bm{k}]\in\mathbb{C}$ is a vector of noise. In this case, we pose recovery as

[TABLE]

where $\delta>0$ is an estimate of the $\ell^{2}$-norm of the error $\|\eta\|$, and $\mathcal{P}_{\Omega}$ denotes projection onto $\Omega$. We make no assumptions on the statistics of the noise $\eta$. In particular, $\eta$ can represent errors due to model-mismatch, such as when the image is not perfectly piecewise constant, or when the edge set of the image does not coincide perfectly with the zero level-set of a bandlimited function.

The following theorem shows that when the deviation of $\hat{f}_{n}$ from $\hat{f}$ is small, the modified recovery program (18) recovers a solution that is close in norm to $\hat{f}$ under the same sampling conditions as Theorem 4.

Theorem 5**.**

Let $f$ be specified by (1), whose edge-set is described by the zero-set of the trigonometric polynomial $\mu_{0}$ bandlimited to $\Lambda_{0}$ with associated incoherence measure $\rho$. Let $\Omega\subset\Gamma$ be an index set drawn uniformly at random within $\Gamma$ such that $|\Omega|$ satisfies the bound (14) in Theorem 4. If the measurements $\hat{f}_{n}$ satisfy $\|\mathcal{P}_{\Omega}({\hat{f}_{n}-\hat{f}})\|_{2}\leq\delta$, then the solution $\hat{g}$ to (18) satisfies

[TABLE]

with probability exceeding $1-|\Gamma|^{-2}$.

See Section IV in the Supplementary Materials for proof. The bound (19) allows us to quantify the effect of model-mismatch on recovery. In particular, suppose the image $f_{n}$ represents a perturbation from an ideal piecewise constant image $f$ such that their difference in $L^{2}$-norm is $\delta$-small:

[TABLE]

Then by Parseval’s theorem, the measurements of $\hat{f}_{n}$ satisfy $\|\mathcal{P}_{\Omega}({\hat{f}_{n}-\hat{f}})\|^{2}\leq\delta$, hence Theorem 5 applies. From (19) we obtain the bound $\|\mathcal{T}(\hat{f}\,)-\mathcal{T}(\hat{g}\,)\|_{F}\leq 5|\Gamma|^{2}\|f_{n}-f\|_{L^{2}}.$ This shows that if the image $f_{n}$ is close to the ideal piecewise constant image $f$ in spatial domain $L^{2}$-norm, then the matrix $\mathcal{T}(\hat{g}\,)$ we recover using (18) will be close in norm to $\mathcal{T}(\hat{f}\,)$ with high probability.

IV Row and column spaces of $\mathcal{T}(\hat{f}\,)$ and incoherence

In this section we define an incoherence measure $\rho$ that satisfies the desired bounds in Prop. 3. We show that the incoherence measure depends only on the geometry of the edge set of the image. The incoherence measure is derived from a new characterization of the row and column spaces of the matrix $\mathcal{T}(\hat{f}\,)$ in terms of Vandermonde-like basis matrices.

IV-A Row and column spaces of $\mathcal{T}(\hat{f}\,)$

Our first lemma gives a basis for the row space of $\mathcal{T}(\hat{f}\,)$:

Lemma 6**.**

A basis of the row space of $\mathcal{T}(\hat{f}\,)$ is given by the columns of the $|\Lambda_{1}|\times R$ Vandermonde-like matrix

[TABLE]

where $\{\bm{k}_{1},...,\bm{k}_{|\Lambda_{1}|}\}$ is a linear indexing of elements in $\Lambda_{1}$, and $P=\{\bm{r}_{1},...,\bm{r}_{R}\}$ is a set of $R=|\Lambda_{1}|-|\Lambda_{1}:\Lambda_{0}|$ distinct points on the edge set curve $\{\mu_{0}=0\}$ chosen such that the columns of $\bm{E}_{\text{row}}$ are linearly independent.

The careful reader will have noticed that Lemma 6 takes for granted the existence of a set of points $P=\{\bm{r}_{1},....,\bm{r}_{R}\}\subset\{\mu_{0}=0\}$ such that the columns of $\bm{E}_{\text{row}}(P)$ is linearly independent. Call such a set $P$ a set of admissible nodes for the curve $\{\mu_{0}=0\}$. The following result shows that sets of admissible nodes always exist and are easy to construct:

Lemma 7**.**

Let $\mu_{0}$ be bandlimited to $\Lambda_{0}$. Any set of $M\geq R+|\Lambda_{0}|$ distinct points on the curve $\{\mu_{0}=0\}$ contains a subset of $R$ points that are a set of admissible nodes.

The next lemma shows that we can characterize the column space of $\mathcal{T}(\hat{f}\,)$ in a similar way as the row space:

Lemma 8**.**

A basis of the column space of $\mathcal{T}(\hat{f}\,)$ is given by the columns of the $2|\Lambda_{2}|\times R$ weighted Vandermonde-like matrix:

[TABLE]

*where where $\{\bm{k}_{1},...,\bm{k}_{|\Lambda_{2}|}\}$ is a linear indexing of elements in $\Lambda_{2}$ and $P=\{\bm{r}_{1},....,\bm{r}_{R}\}$ is a set of admissible nodes for the curve $\{\mu_{0}=0\}$. The weight vectors $\bm{w}_{i}=(w_{i,x},w_{i,y})$, are described by the formula (52) in Appendix VIII, and depend only on the edge set $\{\mu_{0}=0\}$, the nodes $P$, and the filter support $\Lambda_{1}$. *

See Section VIII-C for the proofs of Lemmas 6 and 7, and Section VIII-E for the proof of Lemma 8.

IV-B Incoherence measure

We now show how to define an incoherence measure $\rho$ that satisfies the desired bounds in Prop. 3. Consider the Gram matrix $\bm{G}(P)=[\bm{E}_{\text{row}}(P)]^{*}\bm{E}_{\text{row}}(P)$, where $P$ is any set of $R$ points $\bm{r}_{1},...,\bm{r}_{R}$ on the edge set curve $\{\mu=0\}$. It is easy to see from the definition (21) that the entries of $\bm{G}(P)$ are specified by

[TABLE]

where $D_{\Lambda_{1}}(\bm{r}):=\sum_{\bm{k}\in\Lambda_{1}}e^{j2\pi\bm{k}\cdot\bm{r}}$ is the Dirichlet kernel supported on $\Lambda_{1}$. Note that $\bm{G}(P)$ has ones along the diagonal, and the magnitude of the off-diagonal entries is dictated by the distances $|\bm{r}_{i}-\bm{r}_{j}|$ and the filter support $\Lambda_{1}$. We now define the incoherence measure $\rho$ associated with the edge set $E=\{\mu_{0}=0\}$ in terms of $\bm{G}(P)$.

Definition 9**.**

Suppose the edge set curve $E=\{\mu_{0}=0\}$ has bandwidth $\Lambda_{0}$ (see (3)), and set $R=|\Lambda_{1}|-|\Lambda_{1}\colon\Lambda_{0}|$. Define the incoherence measure $\rho$ by

[TABLE]

where $\lambda_{min}[\bm{G}(P)]$ is the minimum eigenvalue of $\bm{G}(P)$.

Put in words, among all possible arrangements of $R$ points along the edge-set $\{\mu_{0}=0\}$, we seek the arrangement such that the minimum eigenvalue $\bm{G}(P)$ is as large as possible. Intuitively, the optimal arrangement will maximize the minimum separation distance among the $R$ points, and $\rho$ can be thought of as a measure of this geometric property. In particular, edge set curves that enclose a small area, and hence require the points $P$ to be closely spaced along the curve, will result in a large value of $\rho$. According to Theorem 4, the measurement burden will be high for such curves.

Note that curves corresponding to a particular bandwidth can come in different sizes. Specifically, for a fixed $\mu_{0}$ with bandwidth $\Lambda_{0}$ consider the family of curves $\{\mu_{0}=\alpha\}$, where $\alpha$ is a scalar. One can change $\alpha$ to obtain multiple curves with exactly the same bandwidth, each of which correspond to a different levelset of $\mu_{0}$. These level-sets will have different incoherence measures, depending on how large or small the level-set curves are. This shows the incoherence of an edge set captures something besides its bandwidth. See Figure 3 for an illustration.

We can give incoherency measure of an edge set a more precise geometric interpretation based on the minimum separation distance of a set of admissible nodes. We generalize a bound on the condition number of Vandermonde matrices derived in [27] to the case of the Vandermonde-like matrix (21), and use this to derive a bound for the incoherence parameter $\rho$.

Theorem 10**.**

Assume that the points $P=\{(x_{i},y_{i})\}_{i=1}^{R}$ belonging to the curve $\{\mu_{0}=0\}$ satisfy $|x_{i}-x_{j}|>\Delta$ and $|y_{i}-y_{j}|>\Delta$ for all $i\neq j$. Assume the filter support $\Lambda_{1}\subset\mathbb{Z}^{2}$ is a square region symmetric around the origin of size $\sqrt{|\Lambda_{1}|}\times\sqrt{|\Lambda_{1}|}$. Then

[TABLE]

where $\rho$ is the incoherence parameter (24) associated with the curve $\{\mu_{0}=0\}$.

See Section I of the Supplementary Materials for the proof. The bound in (25) shows that the incoherence is close to one (i.e., is as small as possible) when $\Delta\gg 1/\sqrt{|\Lambda_{1}|}$. Since $\Delta$ is the spacing between each pair of points on the curve, to achieve a larger $\Delta$ spacing, and hence a smaller $\rho$, requires a larger curve. This suggests that fewer measurements are required to recover a larger curve, which is consistent with the findings in the isolated Dirac setting [27, 28].

V Numerical Experiments

V-A Algorithms

For small to moderate problem sizes the nuclear norm minimization problem (10) can be solved efficiently with the alternating directions method of multipliers (ADMM) algorithm, which results in a modification of the singular value thresholding (SVT) algorithm [29]. This approach has been proposed for related structured low-rank matrix completion problems in several works, e.g., [30, 3, 11, 12]. We adopt this approach here as well for our small-scale numerical experiments. A detailed implementation of this algorithm can be found in, e.g., [28]. However, we note that for large scale problems, such as those encountered in realistic imaging applications, more efficient approaches need to be adopted, because often in these cases the lifted matrix is too large to be held in memory. A fast algorithm for solving an approximation to (10) for large-scale problems is given in [19].

V-B Phase transitions

In Fig. 4, we study the probability of exact recovery under different assumptions on the filter size and edge set of the image. For these experiments the reconstruction grid $\Gamma$ was of size $65\times 65$. We generated synthetic random piecewise constant functions with known edge set bandwidth (see Fig. 3(c)), and attempted to recover their Fourier coefficients in $\Gamma$ from random samples in $\Omega$ at the specified undersampling factor. For each set of parameters we ran 10 random trials. We count the recovery as “exact” if the recovered coefficients $\hat{f}$ satisfied $\|\hat{f}-\hat{f}_{0}\|/\|\hat{f}_{0}\|<10^{-3}$, where $\hat{f}_{0}$ is the ground truth. The exact recovery rate was then obtained by averaging over the 10 trials.

First, in Fig. 4 (a), we studied the effect of changing the filter size $\Lambda_{1}$ on the recovery while keeping other parameters constant. We fixed the edge-set bandwidth to $|\Lambda_{0}|=9\times 9$ and varied the filter size as $|\Lambda_{1}|=(2K+1)\times(2K+1)$ for $K=1,...,30$. We call $K$ the filter bandwidth. Note that Theorem 4 has restrictions on how large $\Lambda_{1}$ can be. The maximum filter bandwidth for which Theorem 4 holds in this case was $K=15$ (red line in Figure 4(a)), however we extended the filter size to observe the behavior of the algorithm outside of this regime. As predicted by Theorem 4, we find that the optimal performance is obtained when $\Lambda_{1}$ is the largest as allowed by Theorem 4 (roughly half the size of $\Gamma$ in each dimension).

Next, in Fig. 4(b), we study the recovery as a function of the bandwidth of the edge-set of the image. The filter bandwidth was fixed at $K=15$, and we varied the edge-set bandwidth as $|\Lambda_{0}|=(2K_{0}+1)\times(2K_{0}+1)$. The phase transition shows dependence $|\Omega|=O(|\Lambda_{0}|)$ as predicted by Theorem 4.

V-C Comparison with TV minimization on real MRI data

We also compare the proposed Fourier domain interpolation scheme against standard discrete TV minimization in spatial domain:

[TABLE]

Here $\bm{u}\in\mathbb{C}^{N}$ with $N=N_{x}N_{y}$ is a 2-D array representing a discrete $N_{x}\times N_{y}$ image, $\bm{u}_{0}\in\mathbb{C}^{N}$ is the image to be recovered, $\bm{F}\in\mathbb{C}^{N\times N}$ denotes the unitary 2-D discrete Fourier transform (DFT) matrix acting on $N_{1}\times N_{2}$ arrays, $\bm{P}_{\Omega}$ is projection onto the index of sampling locations $\Omega\subset[N_{x}]\times[N_{y}]$, and $TV(\cdot)$ denotes the (isotropic) total variation semi-norm:

[TABLE]

where $\bm{\partial}_{1}$ and $\bm{\partial}_{2}$ are finite difference operators in the horizontal and vertical directions, respectively. The problem (26) has been studied extensively [31, 32, 33, 34, 35, 36] as a model for undersampled MRI reconstruction and other inverse problems in imaging.

In Fig. 5 we perform an experiment comparing against TV minimization and the proposed approach on real MRI data. For this experiment we used a fully-sampled four-coil single-slice acquisition consisting of $256\times 256$ Cartesian k-space samples, which was compressed to a single virtual coil using an SVD-based technique [37]. The data in the single virtual coil was observed to have smoothly varying complex phase in image domain. To compensate for this source of model-mismatch, we further pre-preprocessed the data by removing the complex phase in image domain. We note that this pre-processing step is unrealistic for a true MRI experiment. However, the optimization problem (10) could be modified to incorporate an estimate of the smoothly varying phase in the measurement model; we omit this step for simplicity. Finally, we retrospectively undersampled the pre-processed virtual single coil data, taking $50\%$ uniform random samples. We find that the proposed structured low-rank recovery shows significant improvement recovery error over standard total variation as measured by $\text{SNR}=20\log_{10}(\|\hat{f}\|/\|\hat{f}^{\ast}-\hat{f}\|)$, where $\hat{f}^{*}$ is the recovered data and $\hat{f}$ is the ground truth. The error images indicate the proposed method more faithfully recovers the true edges of the image.

VI Discussion

Discrete domain total-variation minimization has played a central role in compressed sensing from its inception [31, 32], which models the image to be recovered as (approximately) piecewise constant. Since the present work can be thought of as an extension of compressed sensing type guarantees to the continuous domain setting, it is fruitful to explore the connections between our continuous domain model and discrete domain total variation.

At first glance, the structured low-rank matrix completion problem (10) may seem far removed from the TV-minimization problem (26). But, in fact, one can show TV-minimization (26) is equivalent to nuclear norm minimization of a related structured matrix lifting in Fourier domain. Specifically, (26) is equivalent to

[TABLE]

Here

[TABLE]

and $\mathcal{C}_{x}(\bm{F}\bm{u})$, $\mathcal{C}_{y}(\bm{F}\bm{u})$ are block circulant with circulant blocks matrices whose first column is specified by the arrays $\bm{v}_{x}=\bm{F}\bm{\partial}_{x}\bm{u}$ and $\bm{v}_{y}=\bm{F}\bm{\partial}_{x}\bm{u}$. Assuming circular boundary conditions, we can write $(\bm{v}_{x})[k_{x},k_{y}]=(1-e^{j2\pi k_{x}/N_{x}})(\bm{F}\bm{u})[k_{x},k_{y}]$ and $(\bm{v}_{y})[k_{x},k_{y}]=(1-e^{j2\pi k_{y}/N_{y}})(\bm{F}\bm{u})[k_{x},k_{y}]$.

We find it interesting to use this re-formulation of TV-minimization to better understand the proposed approach. In Table I we summarize the similarities and differences. One essential difference is the dimensions of the matrix liftings. In particular, the matrix lifting we propose has dimensions $2|\Lambda_{2}|\times|\Lambda_{1}|$, with $|\Lambda_{1}|\ll|\Lambda_{2}|$ whereas the matrix lifting associated with TV in (28) has dimensions $2N\times N$. If the reconstruction grid size is the same in both cases, i.e., $|\Gamma|=N$, then the proposed matrix lifting has substantially fewer columns than the one associated with TV. This is due to our assumption that edge set of the image has low bandwidth. In other words, we restrict the degrees of freedom of the model by constructing a lifting with fewer columns. We believe this difference may explain the success of the proposed method over TV-minimization observed empirically in Section V.

VII Conclusion

We derived performance guarantees for the recovery of piecewise constant images from random non-uniform Fourier samples via a convex structured low-rank matrix completion problem. This was achieved by adapting results in [3] to the case of a low-rank block two-fold Toeplitz matrix with an additional weighting scheme that arises naturally when considering piecewise constant images. We also define incoherence measures that rely only on geometric properties of the edge set, which indicate that the sampling burden is higher for images with smaller piecewise constant regions.

The recovery guarantees in this work studied the case of uniform random samples. However, in practice we observe that recovery works well with when considering other types of variable density random sampling, where the low spatial frequencies are more heavily sampled. 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.

VIII Appendix A: Incoherence Bounds

VIII-A Notation and Preliminaries

To simplify our arguments, we will convert the linear operators $\mathcal{T}(\hat{f}\,)$ and $\mathcal{T}(\hat{f}\,)^{*}$ defined in Fourier domain to linear operators acting on spaces of trigonometric polynomials (3) in spatial domain. Specifically, for any index set $\Omega\subset\mathbb{Z}^{2}$, let $B_{\Omega}$ denote the vector space of all trigonometric polynomials that have coefficients supported within $\Omega$. Similarly, we denote the space of vector fields $\bm{\rho}=(\rho_{1},\rho_{2})$ with components $\rho_{1},\rho_{2}\in B_{\Omega}$ as $B^{2}_{\Omega}$. We set $\mathcal{S}(f)=\mathcal{F}\mathcal{T}(\hat{f}\,)\mathcal{F}^{-1}$, where $\mathcal{F}$ is the Fourier transform of a periodic function on $[0,1]^{2}$. For any index set $\Lambda$, define the Dirichlet kernel $D_{\Lambda_{1}}(\bm{r}):=\sum_{\bm{k}\in\Lambda_{1}}e^{j2\pi\bm{k}\cdot\bm{r}}$. For all $\varphi\in B_{\Lambda_{1}}$, the action of the linear operator $\mathcal{S}(f):B_{\Lambda_{1}}\rightarrow B_{\Lambda_{2}}^{2}$ can be expressed compactly as

[TABLE]

where $\varphi\,\nabla f$ is understood as a tempered distribution, and the convolution is applied separately to each vector field component. Here convolution with $D_{\Lambda_{2}}$ is a bandlimiting operation. Simliarly, for $\bm{\rho}=(\rho,\rho_{2})\in B^{2}_{\Lambda_{2}}$, the adjoint $\mathcal{S}(f)^{*}$ acts as

[TABLE]

which is the spatial domain equivalent of the adjoint matrix $\mathcal{T}(\hat{f}\,)^{*}$. More expliclty, if $f=1_{U}$ where $U$ is a simply connected region with smooth boundary $\partial U$, a straightforward argument using the divergence theorem shows that the function $\mathcal{S}(f)\varphi$ is given pointwise as the weighted curve integral

[TABLE]

for all $\bm{r}\in[0,1]^{2}$, where $\bm{n}(\bm{r}^{\prime})$ is the outward unit normal to the curve $\partial U$ at $r^{\prime}$, and $ds$ is the arc-length element. Likewise, $\mathcal{S}(f)^{*}\bm{\rho}$ is the function given pointwise by

[TABLE]

for all $\bm{r}\in[0,1]^{2}$. These formulas can be generalized to an arbitrary piecewise constant function $f=\sum_{i}a_{i}1_{U_{i}}$ by linearity. However, in the remainder we focus on the case where $f=1_{U}$ to simplify our arguments.

VIII-B Fundamental subspaces of $\mathcal{S}(f)$ and dimensions

Under the conditions of Theorem 1, the nullspace of $\mathcal{T}(\hat{f}\,)$ is spanned by shifts of the minimal annihilating filter, $\widehat{\mu_{0}}$. In spatial domain, this space consists of all multiples of the minimal degree polynomial $\gamma=\eta\,\mu_{0}$ such that $\gamma$ is bandlimited to $\Lambda_{1}$. We denote this space by

[TABLE]

Note that $(\mu_{0})_{\Lambda_{1}}$ is a subspace of $B_{\Lambda_{1}}$ with dimension $|\Lambda_{1}\colon\Lambda_{0}|$. Therefore, the dimension of the kernel of $\mathcal{S}(f)$, denoted by $\ker\mathcal{S}(f)$, is given by

[TABLE]

By the rank-nullity theorem, the dimension of the image of $\mathcal{S}(f)$, denoted by $\text{im}\,\mathcal{S}(f)$, is

[TABLE]

Likewise, the dimension of the coimage $\text{im}\,\mathcal{S}(f)^{*}$ is also $R$. Furthermore, since $\text{im}\,\mathcal{S}(f)^{*}=[\ker S(f)]^{\perp}$, we have

[TABLE]

This means that any $\gamma\in B_{\Lambda_{1}}$ is in the row space if and only if $\gamma$ is orthogonal to every trigonometric polynomial of the form $\eta\,\mu_{0}\in B_{\Lambda_{1}}$, or equivalently,

[TABLE]

for all $\eta\in B_{\Lambda_{1}\colon\Lambda_{0}}$.

VIII-C Basis for the coimage of $\mathcal{S}(f)$ (corresponding to the row space of $\mathcal{T}(\hat{f})$)

Let $\bm{s}\in[0,1]^{2}$, and set $\varphi_{\bm{s}}\in B_{\Lambda_{1}}$ to be the translated Dirichlet kernel:

[TABLE]

Equivalently, $\varphi_{\bm{s}}\in B_{\Lambda_{1}}$ is the trigonometric polynomial specified in Fourier domain as

[TABLE]

Observe that the inner product of $\varphi_{\bm{s}}$ with any other trigonometric polynomial $\eta\in B_{\Lambda_{1}}$ is given by the point-evaluation of $\eta$ at $\bm{s}$:

[TABLE]

Suppose now that the point $\bm{s}$ satisfies $\mu_{0}(\bm{s})=0$. In this case, we see that $\varphi_{\bm{s}}$ is necessarily in the coimage $\text{im}\,\mathcal{S}(f)^{*}=(\mu_{0})_{\Lambda_{1}}^{\perp}$ since we have

[TABLE]

for any multiple of the minimal polynomial $\gamma\mu_{0}\in B_{\Lambda_{1}}$, i.e., any element in $\ker\mathcal{S}(f)=(\mu_{0})_{\Lambda_{1}}$.

We will now show how to construct a basis for the coimage of $\mathcal{S}(f)$ out of elements having the form $\varphi_{\bm{r}_{i}}$ for some $\bm{r}_{i}$, $i=1,...,R$ belonging to the zero set of $\mu_{0}$. For an arbtirary collection of $R$ points $\{\bm{r}_{i}\}_{i=1}^{R}\subset\{\mu_{0}=0\}$, we are not guaranteed that the set of functions $\{\varphi_{\bm{r}_{i}}\}_{i=1}^{R}$ is linearly independent. However, we will show that there exists a constant $M=M(\Lambda_{0},\Lambda_{1})$ such that for any $M$ distinct points $\{\bm{r}_{i}\}_{i=1}^{M}\subset\{\mu_{0}=0\}$ we can always find a subset of $R$ linearly independent basis functions from the collection $\{\varphi_{\bm{r}_{i}}\}_{i=1}^{M}$. The constant $M$ is related the maximum number of isolated zeros that a system of two trigonometric polynomials can have. The following lemma, which is a consequence of the BKK bound in enumerative algebraic geometry (see, e.g., [38]), puts a bound on $M$. See section II of the supplementary material for proof.

Lemma 11**.**

Let $\Lambda_{1}$ and $\Lambda_{0}$ be rectangular index sets such that $\Lambda_{0}\subset\Lambda_{1}$, and set $R=|\Lambda_{1}|-|\Lambda_{1}\colon\Lambda_{0}|$. For any $\mu_{0},\mu_{1}$ trigonometric polynomials bandlimited to $\Lambda_{0}$ and $\Lambda_{1}$, respectively, the maximum number $M$ of isolated solutions of $\mu_{0}(\bm{r})=\mu_{1}(\bm{r})=0$ is bounded as

[TABLE]

We now prove equivalents of Lemma 6 and Lemma 7 in terms of the spatial domain operator $\mathcal{S}(f)$:

Lemma 12**.**

Let $\{\bm{r}_{1},....,\bm{r}_{N}\}$ be any collection of $N$ distinct points on the curve $\{\mu_{0}=0\}$, where $N\geq R+|\Lambda_{0}|$. Then the coimage space $\text{im}\,\mathcal{S}(f)^{*}=(\mu_{0})^{\perp}_{\Lambda_{1}}$ is spanned by the set of shifted Dirichlet kernels $\varphi_{i}(\bm{r})=D_{\Lambda_{1}}(\bm{r}-\bm{r}_{i})$ for all $i=1,...,N$, i.e.,

[TABLE]

In particular, there exists a subset of $R=|\Lambda_{1}|-|\Lambda_{1}\colon\Lambda_{0}|$ elements of $\left\{\varphi_{\bm{r}_{i}}\right\}_{i=1}^{N}$ that is a basis for the coimage space $(\mu_{0})_{\Lambda_{1}}^{\perp}$.

Proof.

All the functions $\varphi_{\bm{r}_{i}}$ are in $(\mu_{0})_{\Lambda_{1}}^{\perp}$ since we have $\left\langle\varphi_{i},\gamma\mu_{0}\right\rangle=\gamma(\bm{r}_{i})\mu_{0}(\bm{r}_{i})=0$ because each $\bm{r}_{i}$ belong to the zero set of $\mu_{0}$. This implies that

[TABLE]

Our focus is on proving (44) with equality. For this, it is sufficient to show that any vector orthogonal to $\text{span}\{\varphi_{i}\}_{i=1}^{N}$ is in $(\mu_{0})_{\Lambda}$. Assume that there is a vector $\eta(\bm{r})\in B_{\Lambda_{1}}$ that is in the orthogonal complement space of $\text{span}\{\varphi_{i}\}_{i=1}^{N}$. This is only possible if

[TABLE]

Therefore, both $\eta$ and $\mu_{0}$ have $N$ zeros in common. By Lemma 11 this is only possible if $\eta$ contains $\mu_{0}$ as a factor. This implies that all vectors in the orthogonal complement space of $\text{span}\{\varphi_{i}\}_{i=1}^{M}$ are in $(\mu_{0})_{\Lambda_{1}}$, or equivalently

[TABLE]

which together with (45) proves (44).

Finally, we also know that the dimension of $(\mu_{0})_{\Lambda_{1}}^{\perp}$ is equal to $R<M$. Thus, one can select a subset of $R$ basis functions $\varphi_{i}$ that are linearly independent and hence a basis for $(\mu_{0})^{\perp}_{\Lambda_{1}}$.

∎

Translating this result to Fourier domain, we see that the row space of $\mathcal{T}(\hat{f}\,)$ is spanned by the vectors of Fourier coefficients $(\widehat{\varphi_{i}}[\bm{k}]:\bm{k}\in\Lambda_{1})\in\mathbb{C}^{|\Lambda_{1}|}$, for $i=1,...,R$. Equivalently, this can be expressed as the columns of the Vandermonde-like matrix $\bm{E}_{\text{row}}$ specified by (21), which proves Lemma 6 and 7.

VIII-D Discretization of curve integrals: quadrature formula

Using the results from the previous subsection, we now introduce a quadrature formula for curve integrals, which we will use to determine the range space $\text{im}\,\mathcal{S}(f)$ in the next subsection.

Let $\gamma$ be any function in $B_{\Lambda}$ for any $\Lambda\supseteq\Lambda_{0}$. Then from the orthogonal decomposition $B_{\Lambda}=(\mu_{0})_{\Lambda}\oplus(\mu_{0})^{\perp}_{\Lambda}$ we can decompose $\gamma$ as

[TABLE]

where $S=|\Lambda|-|\Lambda\colon\Lambda_{0}|$, and where $\{D_{\Lambda}(\bm{r}-\bm{r}_{i})\}_{i=1}^{S}$ defines a basis of $(\mu_{0})_{\Lambda}^{\perp}$. Here, the coefficients $a_{i}$ in (48) are obtained uniquely as

[TABLE]

where $\bm{D}\in\mathbb{R}^{S\times S}$ is the symmetric matrix with entries $[\bm{D}]_{i,j}=D_{\Lambda}(\bm{r}_{i}-\bm{r}_{j})$ for $1\leq i,j\leq S$. The above expression can be compactly expressed as $\bm{a}=\bm{D}^{-1}\bm{g}$, where $\bm{g}=(\gamma(\bm{r}_{1}),...,\gamma(\bm{r}_{S}))^{T}$.

Lemma 13**.**

Let $f=1_{U}$ where $U$ is a simply connected region with smooth boundary $\partial U$, which is the zero levelset of $\mu_{0}\in B_{\Lambda_{0}}$ and let $\gamma\in B_{\Lambda}$. Consider the curve integral of the form

[TABLE]

where $\bm{n}(\bm{r})=\nabla f(\bm{r})/|\nabla f(\bm{r})|$ is the unit normal on the curve $\partial U$. The curve integral can be evaluated using the quadrature formula

[TABLE]

where the $S=|\Lambda|-|\Lambda:\Lambda_{0}|$ points $\{\bm{r}_{i}\}_{i=1}^{S}$ belong to the curve $\{\mu_{0}=0\}$, and the cooresponding weight vectors $\bm{w}_{i}\in\mathbb{R}^{2}$, $i=1,..,S$, are specified by

[TABLE]

where $\bm{v}_{i}=\oint_{\partial U}D_{\Lambda}(\bm{r}-\bm{r}_{i})\bm{n}(\bm{r})ds(\bm{r})\in\mathbb{R}^{2}$.

Proof.

Decomposing $\gamma(r)$ using (48), we obtain

[TABLE]

The above sum can be expressed in the vector form as

[TABLE]

where $\bm{V}=[\bm{v}_{1}^{T},...,v_{S}^{T}]^{T}\in\mathbb{C}^{R\times 2}$. Setting $\bm{W}=\bm{D}^{-1}\bm{V}=[\bm{w}_{1}^{T},...,\bm{w}_{S}^{T}]^{T}\in\mathbb{C}^{R\times 2}$ we obtain (51). ∎

VIII-E Basis for the range of $\mathcal{S}(f)$ (corresponding to the column space of $\mathcal{T}(\hat{f}\,)$)

We now introduce a basis set for $\mathsf{im}\,\mathcal{S}(f)$, which will be used to prove Lemma 8.

Lemma 14**.**

The range of $\mathcal{S}(f)$, denoted by $\mathsf{im}\,\mathcal{S}(f)$ is specified by

[TABLE]

for an appropriate choice of points $\{\bm{r}_{i}\}_{i=1}^{R}\subset\{\mu_{0}=0\}$ with $R=|\Lambda_{1}|-|\Lambda_{1}\colon\Lambda_{0}|$, and where the weight vectors $\bm{w}_{i}$ are specified by (52).

Proof.

Consider an arbitrary element $\bm{\rho}=(\rho,\rho_{2})\in\mathsf{im}\,\mathcal{S}(f)$. We can express $\bm{\rho}$ as $\bm{\rho}=\mathcal{S}(f)\psi=\mathcal{B}_{\Lambda_{2}}(\psi\nabla f)=D_{\Lambda_{2}}\ast(\psi\nabla f)$ for some $\psi\in B_{\Lambda_{1}}$. By the definition in (33), we have

[TABLE]

where we Lemma 13 in the last step with $S=|\Gamma|-|\Gamma:\Lambda_{0}|$ since the integrand ${\psi(\bm{s})}^{*}\,D_{\Lambda_{2}}(\bm{r}-\bm{s})$ belongs to $B_{\Gamma}$. The above relation shows that any $\bm{\rho}(\bm{s})\in\mathsf{im}\,\mathcal{S}(f)$ can be expressed as the linear combination of the functions $D_{\Lambda_{2}}(\bm{s}-\bm{r}_{i})\,\bm{w}_{i},$ for $i=1,..,S$. Thus, we have $\mathsf{im}\,\mathcal{S}(f)\subset\text{span}\{D_{\Lambda_{2}}(\bm{r}-\bm{r}_{i})\,\bm{w}_{i}\}_{i=1}^{S}$. We also know that $\dim\left(\mathsf{im}\,\mathcal{S}(f)\right)=R<S$. This implies that we can select a subset of $R$ vectors from the set $\{D_{\Lambda_{2}}(\bm{r}-\bm{r}_{i})\,\bm{w}_{i}\}_{i=1}^{S}$ that are linearly independent, which will span $\mathsf{im}\,\mathcal{S}(f)$, and hence define a basis. ∎

Correspondingly, the column space of $\mathcal{T}(\hat{f})$ is spanned by the Fourier coefficients of the basis vectors $\bm{w}_{i}D_{\Lambda_{2}}(\bm{r}-\bm{r}_{i})$, or the columns of the $2|\Lambda_{2}|\times R$ weighted Vandermonde-like matrix $\bm{E}_{\text{col}}$ specified by (22).

VIII-F Incoherence Bounds

VIII-F1 Projection onto row subspace

Let $\bm{E}_{\text{row}}=\bm{E}_{\text{row}}(P)$ be any basis for the row space $V$ of $\mathcal{T}(\hat{f}\,)$ specified by (21), whose columns are vectorized Fourier coefficients of the translated and normalized Dirichlet kernels $\varphi_{i}(\bm{r})=\frac{1}{\sqrt{|\Lambda_{1}|}}D_{\Lambda_{1}}(\bm{r}-\bm{r}_{i})$, $i=1,...,R$, for some set of admissible nodes $P=\{\bm{r}_{1},...,\bm{r}_{R}\}\subset\{\mu_{0}=0\}$. Projecting the measurement basis matrix $\bm{A}_{\bm{k}}$ onto $V$, we have

[TABLE]

Since $\bm{A}_{\bm{k}}$ selects $|\omega(\bm{k})|$ rows of $\bm{E}_{\text{row}}$, each of which has $R$ entries of magnitude $1/\sqrt{|\Lambda_{1}|}$, we have

[TABLE]

where $c_{s}=|\Gamma|/|\Lambda_{1}|$. Hence,

[TABLE]

Minimizing over all sets of admissible nodes $P$ in the construction of $\bm{E}_{\text{row}}$ gives the final bound

[TABLE]

VIII-F2 Projection onto column space

Let $\bm{E}_{\text{col}}=\bm{E}_{\text{col}}(P)$ be a basis for the column space of $\mathcal{T}(\hat{f}\,)$ specified by (22), whose columns are vectorized Fourier coefficients of the translated and weighted Dirichlet kernels $\frac{1}{\sqrt{|\Lambda_{2}|}}\frac{\bm{w}_{i}}{\|\bm{w}_{i}\|}D_{\Lambda_{2}}(\bm{r}-\bm{r}_{i})$, for some set of admissible nodes $P=\{\bm{r}_{1},...,\bm{r}_{R}\}\subset\{\mu_{0}=0\}$. Observe the columns of $\bm{E}_{\text{col}}$ are defined to have unit $\ell^{2}$-norm. Following the same steps as in the row space bound, we have

[TABLE]

Expanding the norm $\|\bm{E}_{\text{col}}^{*}\bm{A}_{\bm{k}}\|_{F}^{2}$ gives

[TABLE]

Hence, we have

[TABLE]

where $\rho^{\prime}$ is defined similarly to $\rho$ as:

[TABLE]

Finally, we show how to bound $\rho^{\prime}$ by $\rho$ in (60). Observe that we can re-define $\rho$ and $\rho^{\prime}$ in terms of the minimum singular value of the basis matrices $\bm{E}_{\text{row}}(P)$ and $\bm{E}_{\text{col}}(P)$, according to the correspondences:

[TABLE]

We will show $\sigma^{2}_{min}(\bm{E}_{\text{row}}(P))\leq\sigma^{2}_{min}(\bm{E}_{\text{col}}(P))$, or equivalently, $[\lambda_{min}(\bm{E}_{\text{col}}(P)^{*}\bm{E}_{\text{col}}(P))]^{-1}\leq[\lambda_{min}(\bm{\bm{E}_{\text{row}}}(P)^{*}\bm{E}_{\text{row}}(P)]^{-1}$, for any set $P$ consisting of $R$ points on the edge set. The claim then follows immediately by taking the minimum over all such sets $P$.

To ease notation, we drop the dependence on the set $P$ in the following. Observe that we can express $\bm{E}_{\text{col}}$ as

[TABLE]

where $\bm{W}_{x}=\text{diag}(\frac{w_{1,x}}{\|\bm{w}_{1}\|},...,\frac{w_{R,x}}{\|\bm{w}_{R}\|})$, $\bm{W}_{y}=\text{diag}(\frac{w_{1,y}}{\|\bm{w}_{1}\|},...,\frac{w_{R,y}}{\|\bm{w}_{R}\|})$, and $\tilde{\bm{E}}_{\text{col}}\in\mathbb{C}^{|\Lambda_{2}|\times R}$ is the Vandermonde-like matrix given entrywise by $[\tilde{\bm{E}}_{\text{col}}]_{i,j}=e^{j2\pi\bm{k}_{i}\cdot\bm{r}_{i}},~{}\text{for all}~{}~{}\bm{k}_{i}\in\Lambda_{2},1\leq j\leq R$. In other words, $\tilde{\bm{E}}_{\text{col}}$ has the same structure as $\bm{E}_{\text{row}}$, but is built with respect to $\Lambda_{2}$ instead of $\Lambda_{1}$. In particular, since we always assume $\Lambda_{1}\subset\Lambda_{2}$, the matrix $\bm{E}_{\text{row}}$ can be embedded as a submatrix of $\tilde{\bm{E}}_{\text{col}}$ by restricting the rows of $\tilde{\bm{E}}_{\text{col}}$ to those indexed by $\Lambda_{1}$. By the variational characterization of the minimum singular value of a matrix, we have

[TABLE]

Finally, since $\bm{E}_{\text{row}}$ is a submatrix of $\tilde{\bm{E}}_{\text{col}}$, we also have $\sigma^{2}_{min}(\bm{E}_{\text{row}})\leq\sigma^{2}_{min}(\tilde{\bm{E}}_{\text{col}})$, which together with (63) gives the desired inequality.

IX Appendix B: Proof of Main Theorem

IX-A Reformulation in lifted domain

We now reformulate the recovery of $\hat{f}$ as a matrix recovery problem in the lifted domain. The matrices $\mathcal{T}_{x}(\hat{f}\,)$ and $\mathcal{T}_{y}(\hat{f}\,)$ contain several copies of the weighted entries $k_{x}\widehat{f}[\bm{k}]$ and $k_{y}\widehat{f}[\bm{k}]$, respectively. We use $\omega(\bm{k})$ to denote the set of locations $(\alpha_{1},\alpha_{2})$ in the matrix $\mathcal{T}_{x}(\hat{f}\,)$ or $\mathcal{T}_{y}(\hat{f}\,)$ that contain the entry $k_{x}\widehat{f}[\bm{k}]$ or $k_{y}\widehat{f}[\bm{k}]$ (this set is the same in either case).

We define the sampling matrices $\bm{A}_{\bm{k}}=\begin{bmatrix}\bm{A}_{1,\bm{k}}\\ \bm{A}_{2,\bm{k}}\end{bmatrix}\in\mathbb{C}^{2|\Lambda_{2}|\times|\Lambda_{1}|}$, for each $\bm{k}=(k_{1},k_{2})\in\Gamma$, where

[TABLE]

for $i=1,2$. The matrices $\{\bm{A}_{\bm{k}}\}_{\bm{k}\in\Gamma}$ form an orthonormal basis for the space of matrices defined by the range of the matrix lifting $\mathcal{T}$; we will call any matrix in the range of $\mathcal{T}$ a structured matrix. For any set of coefficients $\{\widehat{g}[\bm{k}]\}_{\bm{k}\in\Gamma}$ we can expand the structured matrix $\mathcal{T}(\widehat{g})$ as

[TABLE]

We denote the projection operator corresponding to a single sampling location $\bm{k}$ by $\mathcal{A}_{\bm{k}}(\bm{X})=\left\langle\bm{A}_{\bm{k}},\bm{X}\right\rangle\bm{A}_{\bm{k}}$. Since $\{\bm{A}_{\bm{k}}\}_{\bm{k}\in\Gamma}$ is an orthonormal basis, for any structured matrix $\bm{X}$, we have $\sum_{k\in\Gamma}\mathcal{A}_{\bm{k}}(\bm{X})=\mathcal{A}(\bm{X})=\bm{X}$. Since $\bm{A}_{\bm{k}}$ is not the basis for a general $\bm{X}\in\mathbb{C}^{2|\Lambda_{2}|\times|\Lambda_{1}|}$, we also define the projection operator to the space orthogonal to the space of structured matrices by $\mathcal{A}^{\perp}(\bm{X})=(\mathcal{I}-\mathcal{A})(\bm{X})$, where $\mathcal{I}$ is the identity operator. In particular, the constraint $\mathcal{A}^{\perp}(\bm{X})=\mathbf{0}$ implies that $\bm{X}$ is a structured matrix.

The recovery of $f$ from its partial Fourier samples $\widehat{f}[\bm{k}],\bm{k}\in\Omega,$ can thus be reformulated as the completion of a structured matrix $\bm{X}$ from its measurements $\mathcal{A}_{\bm{k}},\bm{k}\in\Omega$. Since the matrix is structured, we have $\mathcal{A}^{\perp}(\bm{X})=\bm{0}$. We thus reformulate (10) as the structured low-rank recovery problem:

[TABLE]

where $\mathcal{Q}_{\Omega}$ that satisifies $\mathbb{E}[\mathcal{Q}_{\Omega}]=\mathcal{I}$ is defined as:

[TABLE]

IX-B Conditions for perfect recovery

The tangent space $T$ of the matrix $\mathbf{X}$ is defined as $T:=\{\mathbf{U}\mathbf{X}_{1}^{H}+\mathbf{X}_{2}\mathbf{V}^{H}:\mathbf{X}_{1}\in\mathbb{C}^{|\Lambda_{2}|\times R},\mathbf{X}_{2}\in\mathbb{C}^{|\Lambda_{1}|\times R}\}$ where $\mathbf{X}=\mathbf{U}\bm{\Lambda}\mathbf{V}^{H}$ is the singular value decomposition of $\bm{X}$. The orthogonal complement of $T$ is denoted by $T^{\perp}$. We first show that if $\mathcal{P}_{T}\approx\mathcal{P}_{T}\mathcal{Q}_{\Omega}\mathcal{P}_{T}$, and if an approximate dual certificate that satisfies certain conditions exist, we obtain perfect recovery.

Lemma 15**.**

Consider a multiset $\Omega$ that contains $m$ random indices. Suppose the sampling operator $\mathcal{Q}_{\Omega}$ obeys

[TABLE]

and there exists a dual certificate matrix $\bm{W}$ satisfying

[TABLE]

Then, $\mathcal{T}(\hat{f})$ is the unique solution to (68), where $n=|\Gamma|$ and $m=|\Omega|$.

See Section III-A of supplementary material for proof. Equation (70) suggests that $\mathcal{Q}_{\Omega}\approx\mathcal{I}$ on the tangent space. The conditions (71), (72), and (73) indicates the existence of a $\bm{W}$, which approximates the exact dual certificate $\bm{U}\bm{V}^{*}$. The above lemma is in line with [3, lemma 1], with the exception of the third condition, indicated by (72). To satisfy (70), we bound the deviation of $\mathcal{P}_{T}\mathcal{Q}_{\Omega}\mathcal{P}_{T}$ from $\mathcal{P}_{T}$ in the following lemma.

Lemma 16**.**

Suppose (12) holds. Then we have

[TABLE]

with probability exceeding $1-n^{-4}$, provided that $m>c_{1}\rho R\,c_{s}\,\log(n)$.

We prove this using [39, Theorem 1.6]. (See Section III-B of supplementary material)

IX-C Construction of the approximate dual certificate $\bm{W}$

We will now use the golfing scheme of [25, 3] to construct an approximate dual certificate $\bm{W}$, which satisfies (71), (72), and (73). In particular, we generate $j_{0}$ independent random sampling sets $\Omega_{i};1\leq i\leq j_{0}$, each containing $\tilde{m}=m/j_{0}$ samples corresponding to sampling with replacement. We start with $\bm{F}_{0}=\bm{U}\bm{V}^{*}$, and follow the following steps:

$~{}\bm{F}_{0}=\bm{U}\bm{V}^{*}$ and set $j_{0}=3\log_{\frac{1}{\epsilon}}n$. 2. 2.

$~{}\forall i(1\leq i\leq j_{0}),\bm{F}_{i}=\mathcal{P_{T}}(\mathcal{I}-\mathcal{Q}_{\Omega_{i}})\mathcal{P}_{T}(\bm{F}_{i-1})$ 3. 3.

$~{}\bm{W}=\sum_{j=1}^{j_{0}}\mathcal{Q}_{\Omega_{i}}\bm{F}_{j-1}$

Step 3 ensures that $\bm{W}$ satisfies (71) since each term $\bm{W}_{i}=\mathcal{Q}_{\Omega_{i}}\bm{F}_{j-1}$ satisfies $\mathcal{Q}_{\Omega}^{\perp}(\bm{W}_{i})=0$. The recursive construction also satisfies (72). In particular,

[TABLE]

Now we focus on showing that $\bm{W}$ satisfies (73). Note that if $j_{0}$ is chosen as $3\log_{\frac{1}{\epsilon}}n$, assuming $n>6$, we have $\left(\epsilon\right)^{j_{0}}\;n<\frac{1}{6~{}n}$.

Lemma 17**.**

For any matrix $\bm{M}$, there exists some numerical constant $c_{2}$ such that

[TABLE]

with probability at least $1-n^{-10}$. Here,

[TABLE]

See Section III-C of supplementary material for proof.

Lemma 18**.**

Assume that there exists a constant $\mu_{5}$ such that $\omega_{\bm{k}}\|\mathcal{P}_{T}(\bm{A}_{\bm{k}})\|_{\mathcal{A},2}\leq\frac{\mu_{5}R}{n}$. For any matrix $\bm{M}$, we have

[TABLE]

with probability at least $1-n^{-10}$.

See Section III-D of the Supplementary Materials for proof.

Lemma 19**.**

For any matrix $\bm{M}\in T$, there exists some numerical constant $c_{4}$, such that

[TABLE]

with probability at least $1-n^{-10}$.

See Section III-F of the Supplementary Materials for proof.

From the golfing scheme, we have $\|\mathcal{P}_{T^{\perp}}(\bm{W})\|\leq\sum_{j=1}^{j_{0}}\|\mathcal{P}_{T^{\perp}}\mathcal{Q}_{\Omega_{i}}\mathcal{P}_{T}\bm{F}_{j-1}\|$. Using lemma 17 and substituting from lemma 18 and lemma 19, we have

[TABLE]

The last inequality holds if $\tilde{m}=m/j_{0}\gg\max\left(\mu_{5},\rho c_{s}\right)R\log n$. Substituting for $j_{0}=3\log_{\frac{1}{\epsilon}}(n)$ assumed in the golfing scheme, we require $m\gg c_{6}~{}\max\left(\mu_{5},\rho c_{s}\right)R\log^{2}n$ to satisfy the above inquality. See Section III-G of the Supplementary Materials for details. We will now present the lemmas bounding $\|\bm{F}_{0}\|_{\mathcal{A},2}$ and $\|\bm{F}_{0}\|_{\mathcal{A},\infty}$, where $\bm{F}_{0}=\bm{U}\bm{V}^{*}$.

Lemma 20**.**

With the incoherence measure $\rho$, one can bound

[TABLE]

for $\mu_{3}=3\rho$ and $c_{7}$ is some constant.

See Section III-H of the Supplementary Materials for proof. From (83), we see that the constant $\mu_{5}$ in lemma 19 can be chosen as $\mu_{5}=c_{7}~{}\mu_{3}~{}c_{s}~{}\log^{2}(n)$ such that $\omega_{\bm{k}}\|\mathcal{P}_{T}(\bm{A}_{\bm{k}})\|_{\mathcal{A},2}\leq\frac{\mu_{5}R}{n}$. Substituting for $\mu_{5}$, we observe that the dominant term has its dependence on $\log^{4}(n)$. Thus, $\|\mathcal{P}_{T^{\perp}}\mathcal{Q}_{\Omega_{i}}\bm{F}_{j-1}\|<1/2$ if $m>c_{6}c_{7}\,c_{s}\,\left(3\rho\right)\,R\,\log^{4}(n).$

I Proof of Theorem 10

The proof we give is multi-dimensional generalization of the proof of [27]. We will make use the following lemma from [27]:

Lemma 21**.**

There is an entire function $c_{E}(t)$ whose Fourier transform is supported in the interval $-\Delta,\Delta$, which satisfies $c_{E}(t)<I_{E}(t)$—the indicator function of the interval $E=[-n/2,n/2]$.

[TABLE]

The above function $c_{E}(t)$ is known as the Beurling-Selberg minorant of $I_{E}(t)$. Note that $\widehat{I_{E}}(0)=n$ and hence $\widehat{c_{E}}(0)=n-1/\Delta$, where $\widehat{I_{E}}$ is the Fourier transform of $I_{E}$.

We now give the proof of Theorem 10:

Proof.

We note that $\lambda_{min}\left(\bm{E}_{\text{row}}^{*}\bm{E}_{\text{row}}\right)=\min_{\|\bm{u}\|=1}\|\bm{E}_{\text{row}}\bm{u}\|^{2}$. From the definition of $\bm{E}_{\text{row}}$, we have

$\bm{E}_{\text{row}}\bm{u}=\frac{1}{\sqrt{|\Lambda_{1}|}}\sum_{i=1}^{R}u_{i}\exp\left(j2\pi\bm{k}\cdot\bm{r}_{i}\right);\bm{k}\in\Lambda_{1}$. We consider the continuous domain function

[TABLE]

and rewrite the discrete summation in $\|\bm{E}_{\text{row}}\bm{u}\|^{2}$ as the integral

[TABLE]

where $h(f)=\sum_{m=-\infty}^{\infty}\exp\left(j2\pi fm\right)$ is the Dirac comb function and $I_{E}(f)$ is the indicator function of the region

$[-\sqrt{\frac{|\Lambda_{1}|}{2}},\sqrt{\frac{|\Lambda_{1}|}{2}}]$. Minorizing $I_{E}$ by $c_{E}$, specified by Lemma 21, we obtain

[TABLE]

In the last step, we used the expression of the Dirac comb function. Since the Fourier transform of $c_{E}$ is supported within $[-\Delta,\Delta]$, the terms $\widehat{c_{E}}(x_{i}-x_{j}-m)$ and $\widehat{c_{E}}(x_{i}-x_{j}-m)$ are non-zero only if $x_{i}=x_{j};m=0$ and $y_{i}=y_{j};m=0$, since $|x_{i}-x_{j}|>\Delta,|y_{i}-y_{j}|>\Delta;i\neq j$. Thus, we have

[TABLE]

∎

II The BKK bound (Proof of Lemma 11)

The BKK bound in enumerative algebraic geometry (see, e.g., [38]) is a well-known result that relates maximum number of isolated solutions of a system of polynomials to their coefficient support sets. Specifically, the BKK bound is typically stated in terms of Laurent polynomials, i.e., functions of the form $q(t,s)=\sum_{k=-\infty}^{\infty}\sum_{\ell=-\infty}^{\infty}t^{k}s^{\ell}$ with $t,s\in\mathbb{C}$, and only finitely many non-zero coefficients $c_{k,\ell}\in\mathbb{C}$. Since a trigonometric polynomial is the restriction of a Laurent polynomial to the complex unit torus $\{(t,s)=(e^{j2\pi x},e^{j2\pi y}):x,y\in[0,1)\}$, the result also holds for trigonometric polynomials, which we state below:

Theorem** (BKK Bound).**

Let $\mu_{1}$ and $\mu_{2}$ be trigonometric polynomials with coefficient supports $\Omega_{1}$ and $\Omega_{2}$, and let $P_{1}=\text{conv}(\Omega_{1})$ and $P_{2}=\text{conv}(\Omega_{1})$, where $\text{conv}(\cdot)$ denotes the convex hull of a set in $\mathbb{Z}^{2}$ treated as a subset of $\mathbb{R}^{2}$. The number of isolated solutions of the system $\mu_{1}(\bm{r})=\mu_{2}(\bm{r})=0$ for $\bm{r}\in[0,1)^{2}$ is at most

[TABLE]

where $area(\cdot)$ denotes the usual Euclidean area. In particular, if $\mu_{1}$ is irreducible, and the common isolated zeros of $\mu_{2}$ and $\mu_{1}$ are greater than $\mathcal{M}(\Omega_{1},\Omega_{2})$, then $\mu_{1}$ must divide $\mu_{2}$.

When $\Lambda_{0}$ and $\Lambda_{1}$ are rectangular index sets satisfying $\Lambda_{0}\subset\Lambda_{1}$, a straightforward computation reveals that

[TABLE]

which establishes the bound in Lemma 11.

III Proofs of results in Appendix B

III-A Proof of Lemma 15

Proof.

Let $\bm{X}$ be the unique minimizer for the convex optimization problem (68) and $\bm{H}$ be a perturbation of $\bm{X}$. To prove exact recovery of (68), it suffices to show the existence of an exact dual certificate $\bm{W}$, s.t, $\|\bm{X}+\bm{H}\|_{*}>\|\bm{X}\|_{*}+\left\langle\bm{W},\bm{H}\right\rangle$. With an approach similar to [25], we now show that the existence of an approximate dual certificate will guarantee unique recovery. We separately consider the two cases based on the relative energies of $\bm{H}_{T}$ and $\bm{H}_{T}^{\perp}$, where $\bm{H}=\bm{H}_{T}+\bm{H}_{T}^{\perp}$ and $T$ denotes the tangent space.

III-A1 Case 1: $\|\bm{H}_{T}\|_{F}>2n\,\|\bm{H}_{T}^{\perp}\|_{F}$

We will show that $\bm{X}+\bm{H}$ is infeasible (i.e, $\|\mathcal{Q}_{\Omega}\bm{H}\|_{F}>0$), if $\mathcal{Q}_{\Omega}\approx\mathcal{I}$ on the tangent space (i.e, (70) is satisfied). We have

[TABLE]

We upper bound the second term as $\|\mathcal{Q}_{\Omega}\bm{H}_{T}^{\perp}\|_{F}\leq\|\mathcal{Q}_{\Omega}\|~{}\|\bm{H}_{T}^{\perp}\|_{F}$. By definition, we have $\|\mathcal{Q}_{\Omega}\|=\|~{}\frac{n}{m}~{}\mathcal{A}_{\Omega}+\mathcal{A}^{\perp}\|$, from which we obtain $\left\|\mathcal{Q}_{\Omega}\right\|=\left\|\frac{n}{m}\mathcal{A}_{\Omega}+\mathcal{A}^{\perp}\right\|\leq\|\frac{n}{m}\left(\mathcal{A}+\mathcal{A}^{\perp}\right)\|=\frac{n}{m}<n$. We omit the m term in the denominator to remove the dependence on m. Hence, we have $\|\mathcal{Q}_{\Omega}\bm{H}_{T}^{\perp}\|_{F}\leq n\|\bm{H}_{T}^{\perp}\|_{F}$. We now lower bound the first term in (92):

[TABLE]

Since we assumed that $\|\bm{H}_{T}\|_{F}>n\,\|\bm{H}_{T}^{\perp}\|_{F}$, we have $\|\mathcal{Q}_{\Omega}\bm{H}\|_{F}\geq\left(\sqrt{2}-1\right)\,n~{}\|\bm{H}_{T}^{\perp}\|_{F}>0$, implying that such an $\bm{H}$ is infeasible.

III-A2 Case 2: $\|\bm{H}_{T}\|_{F}\leq 2n\,\|\bm{H}_{T}^{\perp}\|_{F}$

We now show that if $\bm{H}_{T}$ is small and $\bm{X}+\bm{H}$ is feasible, then the nuclear norm of $\bm{X}+\bm{H}$ is larger than $\bm{X}$. Since $\bm{X}+\bm{H}$ is feasible, we have $\mathcal{Q}_{\Omega}(\bm{H})=0.$ If $\bm{X}=\bm{U}\bm{\Sigma}\bm{V}^{*}$ represents the singular value decomposition of $\bm{X}$, the subgradient of $\|\bm{X}\|_{*}$ is parametrized as $\bm{U}\bm{V}^{*}+\bm{Z}_{0};~{}\bm{Z}_{0}\in T^{\perp}$. By definition of subgradient, we have $\|\bm{X}+\bm{H}\|_{*}\geq\|\bm{X}\|_{*}+\left\langle\bm{U}\bm{V}^{*},\bm{H}\right\rangle+\left\langle\bm{Z}_{0},\bm{H}\right\rangle$. We consider a $\bm{W}$ in the range space of $\mathcal{Q}_{\Omega}$ that satisfies (71) $\mathcal{Q}_{\Omega}^{\perp}(\bm{W})=0$, where $\mathcal{Q}_{\Omega}^{\perp}$ is the projection $\mathcal{Q}_{\Omega}^{\perp}=\mathcal{A}-\mathcal{A}_{\Omega}$. Using this $\bm{W}$, we rewrite the above relation as

[TABLE]

The second term in (93) vanishes since $\bm{W}$ lives in the range of $\mathcal{Q}_{\Omega}$ (i.e., $\mathcal{Q}_{\Omega}^{\perp}(\bm{W})=0$) while $\bm{H}$ lives in the kernel of $\mathcal{Q}_{\Omega}$ (i.e., since $\mathcal{Q}_{\Omega}(\bm{H})=0$). Since $\bm{Z}_{0}\in T^{\perp}$ and $\|\bm{Z}_{0}\|\leq 1$, the third term is less than or equal to $\|P_{T}^{\perp}\bm{H}\|_{*}$. We now focus on the last term. If $\bm{W}$ satisfies (72) and (73), we have

[TABLE]

Substituting in (93) and using $\|\bm{M}\|_{*}\geq\|\bm{M}\|_{F}$, we have

[TABLE]

We used $\|\bm{H}_{T}\|_{F}\leq 2n\,\|\bm{H}_{T}^{\perp}\|_{F}$ in the above. The above inequality implies that $\|\bm{X}+\bm{H}\|_{*}\geq\|\bm{X}\|_{*}$ and hence $\bm{X}$ is the unique minimizer. ∎

III-B *Proof of Lemma 16 *

Proof.

Substituting the definition of $\mathcal{Q}_{\Omega}$ in (70), we see that

[TABLE]

We set $\mathcal{Z}_{\bm{k}}=\frac{n}{m}~{}\mathcal{P}_{T}\mathcal{A}_{\bm{k}}\mathcal{P}_{T}$ with $\mathbb{E}[\mathcal{Z}_{\bm{k}}]=\frac{1}{m}~{}\mathcal{P}_{T}\mathcal{A}\mathcal{P}_{T}$, we obtain

[TABLE]

We apply the operator Bernstein’s inequality [39, Theorem 1.6] to determine $\mathbb{P}\left(\left\|\sum_{\bm{k}\in\Omega}\mathcal{S}_{\bm{k}}\right\|\geq\epsilon\right)$. Using steps similar to [3, Lemma 3], we obtain the upper bounds $\left\|\mathbb{E}\left[\mathcal{S}_{k}^{2}\right]~{}\right\|\leq\frac{4\rho Rc_{s}}{m}$ and $\|\mathcal{S}_{k}\|\leq\frac{4\rho Rc_{s}}{m}$. We now apply the matrix Bernstein’s inequality [39, Theorem 1.6] to obtain

[TABLE]

We desire $\mathbb{P}\left(\left\|\sum_{\bm{k}\in\Omega}\mathcal{S}_{\bm{k}}\right\|\geq\epsilon\right)<(n)^{-b}$. Setting these values in the above inequality, we obtain

$(n)^{-b}\geq n~{}\exp\left(\frac{-3m}{112~{}\rho Rc_{s}}\right)$. Taking log of both sides and simplifying, we obtain

[TABLE]

∎

III-C Proof of Lemma 17

Proof.

We define

[TABLE]

which satisfies $\mathbb{E}(\bm{S}_{k})=0$ and $\left\|\mathcal{I}-\mathcal{Q}_{\Omega}\right\|=\|\sum_{\bm{k}\in\Omega}\mathcal{S}_{k}\|$. To bound the right hand side using operator Bernstein’s inequality, we require the bounds $\|\mathcal{S}_{k}(\bm{M})\|\leq B$ and $\sigma^{2}=\max\left\{\left\|\sum_{k}~{}\mathbb{E}[\mathcal{S}_{k}\mathcal{S}_{k}^{*}]\right\|,\left\|\sum_{k}~{}\mathbb{E}[\mathcal{S}_{k}^{*}\mathcal{S}_{k}]\right\|\right\}$. We first consider

[TABLE]

which gives $\left\|\mathbb{E}\left[\sum_{k\in\Omega}~{}\mathcal{S}_{k}^{*}\mathcal{S}_{k}\right]\right\|=\frac{m}{n}\left\|\sum_{k\in\Gamma}~{}\mathcal{S}_{k}^{*}\mathcal{S}_{k}\right\|\leq\frac{n}{m}\|\bm{M}\|_{\mathcal{A},2}^{2}$. Similarly, we have

[TABLE]

where $\bm{k}=(k_{1},k_{2})^{T}$, which also gives the bound $\left\|\mathbb{E}\left[\sum_{k\in\Omega}~{}\mathcal{S}_{k}\mathcal{S}_{k}^{*}\right]\right\|\leq\frac{n}{m}\|\bm{M}\|_{\mathcal{A},2}^{2}$. Here, $\bigotimes$ denotes the Kroneker product. Similar to the arguments in [3], we have $\|\mathcal{S}_{k}\bm{M}\|\leq\frac{2n}{m}\|\bm{M}\|_{\mathcal{A},\infty}$. Combining these terms into [3, Lemma 11], the result is proved. ∎

III-D Proof of Lemma 18

Proof.

We note that

[TABLE]

We assume that $\mathcal{A}_{\Omega}=\sum_{i=1}^{m}\mathcal{A}_{{\bm{\alpha}}_{i}}$, where ${\bm{\alpha}}_{i}$ are independent indices picked at random. Correspondingly, we consider vectors $\bm{z}_{{\bm{\alpha}}_{i}}$ of length $n=|\Gamma|$, whose entries are specified by

[TABLE]

Note that the desired bound $\left\|\mathcal{P}_{T}[(\mathcal{I}-\mathcal{Q}_{\Omega})(\bm{M})]\right\|_{\mathcal{A},2}=\left\|\sum_{i=1}^{m}\bm{z}_{{\bm{\alpha}}_{i}}\right\|_{2}$. We have $\mathbb{E}(\bm{z}_{{\bm{\alpha}}})=0$. We proceed as [3, Lemma 5] and from the definition of $\mu_{5}$, we have, $\|\bm{z}_{{\bm{\alpha}}_{i}}\|_{2}\leq 2\sqrt{\frac{n\mu_{5}R}{m^{2}}}\|\bm{M}\|_{\mathcal{A},\infty}$ and as $\bm{z}_{{\bm{\alpha}}_{i}}$’s are vectors,

[TABLE]

Substituting these bounds in the operator Bernstein’s inequality[3, Lemma 11], the result is proved. ∎

III-E Incoherence between two sampling bases

The proof of Lemma 19 relies on the bound on $\left|\left\langle\bm{A}_{\bm{\beta}},\mathcal{P}_{T}\bm{A}_{\bm{\alpha}}\right\rangle\right|$. We introduce the following lemma to establish the incoherence between two sampling bases:

Lemma 22**.**

Under the incoherence conditions of proposition 3 and definition 9 of the main text,

[TABLE]

We will now bound $\left|\left\langle\bm{A}_{\beta},\mathcal{P}_{T}\bm{A}_{{\bm{\alpha}}}\right\rangle\right|$. We see that

[TABLE]

We will now bound each of the terms in the right hand side. We observe that $\bm{A}_{\beta}$ has $\omega_{\beta}$ entries of magnitude $\frac{k_{\beta,1}}{\|\bm{k}_{\beta}\|\sqrt{2\omega_{\beta}}}$ and $\omega_{\beta}$ entries of magnitude $\frac{k_{\beta,2}}{\|\bm{k}_{\beta}\|\sqrt{2\omega_{\beta}}}$. Hence, we have $\|\bm{A}_{\beta}\|_{\ell_{1}}\leq\sqrt{\omega_{\beta}}$. We consider

[TABLE]

We now bound the entries of $\bm{U}\bm{U}^{*}$:

[TABLE]

which gives $\left|\left\langle\bm{A}_{\beta},\bm{U}\bm{U}^{*}\bm{A}_{{\bm{\alpha}}}\right\rangle\right|\leq\sqrt{\frac{\omega_{\beta}}{\omega_{{\bm{\alpha}}}}}\frac{\rho\,c_{s}\,R}{n}.$ Proceeding along the same lines, we obtain $\left|\left\langle\bm{A}_{\beta},\bm{A}_{{\bm{\alpha}}}\bm{V}\bm{V}^{*}\right\rangle\right|\leq\sqrt{\frac{\omega_{\beta}}{\omega_{{\bm{\alpha}}}}}\frac{\rho\,c_{s}\,R}{n}$ and $\left|\left\langle\bm{U}\bm{U}^{*}\bm{A}_{{\bm{\alpha}}}\bm{V}\bm{V}^{*},\bm{A}_{{\bm{\alpha}}}\right\rangle\right|\leq\sqrt{\frac{\omega_{\beta}}{\omega_{{\bm{\alpha}}}}}\frac{\rho\,c_{s}\,R}{n}$. Substituting these expressions into (101), we prove the Lemma.

III-F Proof of Lemma 19

Proof.

Proceeding with the definition of $\bm{z}_{{\bm{\alpha}}}$ in (100), we observe that

[TABLE]

Using Lemma 22, we obtain

[TABLE]

Substituting (22) in (105), we have $\left|\bm{z}_{{\bm{\alpha}}_{i}}^{k}\right|\leq\frac{2c_{s}r\rho}{m}\|\bm{M}\|_{\mathcal{A},\infty}$. Similarly, we have

[TABLE]

Substituting (105) in (107), we have $|\mathbb{E}[\sum_{i=1}^{m}|\bm{z}^{k}_{{\bm{\alpha}}_{i}}|^{2}]|=\frac{(2c_{s}r\rho)^{2}}{mn}\|\bm{M}\|^{2}_{\mathcal{A},2}$. We have the necessary terms to bound $\left|\sum_{i=1}^{m}\bm{z}_{{\bm{\alpha}}_{i}}^{k}\right|$ for any $k$, which can bound $\max_{\begin{subarray}{c}k\end{subarray}}\left|\sum_{i=1}^{m}\bm{z}_{{\bm{\alpha}}_{i}}^{k}\right|$. We apply [3, Lemma 11], to prove the result. ∎

III-G Upper bound for $\|\mathcal{P}_{T}^{\perp}(\bm{W})\|$

From the golfing scheme, we have $\|\mathcal{P}_{T^{\perp}}(\bm{W})\|\leq\sum_{j=1}^{j_{0}}\|\mathcal{P}_{T^{\perp}}\mathcal{Q}_{\Omega_{i}}\mathcal{P}_{T}\bm{F}_{j-1}\|$. Each of the terms in the right hand side can be bounded as

[TABLE]

We used Lemma 17 in (III-G). Here, $\tilde{m}=m/j_{0}=|\Omega_{j}|$. Substituting from Lemma 18 and Lemma 19, we get

[TABLE]

The last inequality holds if $\tilde{m}=m/j_{0}\gg\max\left(\mu_{5},\rho c_{s}\right)R\log n$. Substituting for $j_{0}=3\log_{\frac{1}{\epsilon}}(n)$ assumed in the golfing scheme, we require

[TABLE]

to satisfy the above inquality. We will now focus on bounding $\|\bm{F}_{0}\|_{\mathcal{A},2}$ and $\|\bm{F}_{0}\|_{\mathcal{A},\infty}$, where $\bm{F}_{0}=\bm{U}\bm{V}^{*}$.

III-H Proof of Lemma 20

Proof.

The proof of this theorem is in line with [3, Lemma 7]. The first term is upper bounded by the maximum entry of the matrix (i.e, $\|\bm{U}\bm{V}^{*}\|_{\mathcal{A},\infty}=\max_{k}|(\bm{U}\bm{V}^{*})_{k}|$).

[TABLE]

We now show that the energy of each row of the matrices $\bm{U}\bm{V}^{*}$ and $\mathcal{P}_{T}\left(\sqrt{\omega_{{\bm{\alpha}}}}~{}\bm{A}_{{\bm{\alpha}}}\right)$ are upper bounded; we use this relation to prove (82) and (83). The energy of the rows of $\bm{U}\bm{V}^{*}$ are bounded as

[TABLE]

Similarly, by the definition of $\mathcal{P}_{T}$, we have

[TABLE]

We now bound each of the terms in the right hand side of the above expression as

[TABLE]

Here, we use the property that the operator norm of $\bm{U}$ and $\sqrt{\omega_{\bm{\alpha}}}\bm{A}_{{\bm{\alpha}}}$ are bounded by 1. Substituting in (112), we obtain the upper bound for the energy of the rows as

[TABLE]

Now, applying [3, Lemma 12] that relates the upper bound of $\|\bm{M}\|_{\mathcal{A},2}^{2}$ to the upper bound of the energy of the rows $\max_{i}\|\bm{e}_{i}\bm{M}\|^{2}$, we obtain the results (82) and (83).

∎

Substituting (81) and (82) in (109), we get

[TABLE]

Similar to the argument before, if $\tilde{m}=m/j_{0}\gg c_{7}\mu_{3}c_{s}R\log^{3}n$ and $\tilde{m}=m/j_{0}\gg\rho c_{s}R\log n$, or equivalently $m\gg{\rho c_{s}}r\log^{2}n$, we have $\|\mathcal{P}_{T^{\perp}}\mathcal{Q}_{\Omega_{i}}\bm{F}_{j-1}\|\leq\left(\frac{1}{2}\right)^{j_{0}}$. Combining these conditions with (110), we need

[TABLE]

From (83), we see that the constant $\mu_{5}$ in Lemma 19 can be chosen as

[TABLE]

such that $\omega_{\bm{k}}\|\mathcal{P}_{T}(\bm{A}_{\bm{k}})\|_{\mathcal{A},2}\leq\frac{\mu_{5}R}{n}$. Substituting for $\mu_{5}$ from (115), we observe that the dominant term is the first one due to its dependence on $\log^{4}(n)$. Thus, $\|\mathcal{P}_{T^{\perp}}\mathcal{Q}_{\Omega_{i}}\bm{F}_{j-1}\|<1/2$ if

[TABLE]

IV Proof of Theorem 5

The proof is similar to that in Appendix L of [3]. The optimization problem (18) can be reformulated in the lifted matrix domain as:

[TABLE]

Let $\widehat{\bm{X}}=\mathcal{T}(\hat{f}\,)+\bm{H}$ be the solution to (117), where $\bm{H}=\mathcal{A}_{\Omega}(\bm{H})+\mathcal{A}_{\Omega^{\perp}}(\bm{H})$. Using (117),

[TABLE]

From the argument in (92), we only consider the case,

[TABLE]

defined in III-A2 Case 2. Now, $\|\mathcal{T}(\hat{f}\,)\|_{*}\geq\|\widehat{\bm{X}}\|_{*}\geq\|\mathcal{T}(\hat{f}\,)+\mathcal{A}_{\Omega^{\perp}}(\bm{H})\|_{*}-\|\mathcal{A}_{\Omega}(\bm{H})\|_{*}$. From (94) and the treatement in Appendix L, [3], we have $\|\mathcal{T}(\hat{f}\,)\|_{*}\geq\|\mathcal{T}(\hat{f}\,)\|_{*}+\frac{1}{6}\|\mathcal{P}_{T}\mathcal{A}_{\Omega^{\perp}}(\bm{H})\|_{F}-\|\mathcal{A}_{\Omega}(\bm{H})\|_{*}$, which gives $\|\mathcal{P}_{T}\mathcal{A}_{\Omega^{\perp}}(\bm{H})\|_{F}\leq 6\|\mathcal{A}_{\Omega}(\bm{H})\|_{*}\leq 6\sqrt{n}\|\mathcal{A}_{\Omega}(\bm{H})\|_{F}\leq 12n\delta$. Finally using (118),

[TABLE]

which concludes the proof.