Image Restoration: Structured Low Rank Matrix Framework for Piecewise Smooth Functions and Beyond

Recently, mapping a signal/image into a low rank Hankel/Toeplitz matrix has become an emerging alternative to the traditional sparse regularization, due to its ability to alleviate the basis mismatch between the true support in the continuous domain and the discrete grid. In this paper, we introduce a novel structured low rank matrix framework to restore piecewise smooth functions. Inspired by the total generalized variation to use sparse higher order derivatives, we derive that the Fourier samples of higher order derivatives satisfy an annihilation relation, resulting in a low rank multi-fold Hankel matrix. We further observe that the SVD of a low rank Hankel matrix corresponds to a tight wavelet frame system which can represent the image with sparse coefficients. Based on this observation, we also propose a wavelet frame analysis approach based continuous domain regularization model for the piecewise smooth image restoration. Finally, numerical results on image restoration tasks are presented as a proof-of-concept study to demonstrate that the proposed approach is compared favorably against several popular discrete regularization approaches and structured low rank matrix approaches.


Introduction
Image restoration, including image denoising, deblurring, inpainting, computed tomography, etc., is one of the central problems in imaging science. It aims at recovering an image of high-quality from a given measurement which is degraded during the process of imaging, acquisition, and communication. The image restoration problem is typically modeled as the following linear inverse problem: where f is the degraded measurement or the observed image, ζ is a certain additive noise, and A is some linear operator which takes different forms for different image restoration problems.
Since the operator A is in general ill-conditioned or non-invertible, it is in general to use a regularization on the images to be restored. Most widely used examples include the variational approaches including total variation (TV) [61] and its nonlocal variants [68], the inf-convolution [20], the total generalized variation (TGV) [6,7], and the combined first and second order total variation [1,55]. Apart from the variational approaches, the applied harmonic analysis approaches including curvelets [17], Gabor frames [42], shearlets [45], complex tight framelets [38], and wavelet frames [8,9,11,21] are also widely used in the literature. The common concept of these methods is to use a sparse regularization on discrete images under a discrete linear transformation to regularize smooth image components while preserving image singularities such as edges, ridges, and corners. However, in many applications, the true singularities lie in the continuous domain, and the discretization error will lead to the basis mismatch [23,52,67] between the true singularities and the discrete image grid. Such a basis mismatch would destroy the sparse structure of the image, and thus can degrade the restoration quality [67].
Recently, continuous domain regularization is emerging as a powerful alternative to the discrete domain sparse regularization [2,19,22,49]. By such an "off-the-grid" approach, we can exploit the sparsity prior in continuous domain, which enables us to alleviate the basis mismatch due to the discretization [49]. To the best of our knowledge, this off-the-grid regularization stems from the Prony's method [58] which corresponds a superposition of a few sinusoids to a structured low rank matrix for the Dirac stream retrieval. Hence, we can adopt the so-called structured low rank matrix (SLRM) approach [22,53,66] for the restoration of spectrally sparse signal whose Fourier transform is a Dirac stream [14,15,16]. Apart from the spectrally sparse signal restoration, the SLRM can also be used to restore Fourier samples of a one dimensional piecewise constant signal [3,63] as in this case the Fourier samples of a derivative becomes the superposition of a few sinusoids. However, even though the SLRM can be easily applied to the case of isolated singularities [18,65], the extension to the multi dimensional (piecewise smooth) image restoration is not straightforward [52]. Since the image singularities such as the edges and ridges in general form a continuous curve in a two dimensional domain, it is in general challenging to construct a structured low rank matrix from the Fourier samples.
In this paper, we introduce a new structured low rank matrix framework for the piecewise smooth image restoration. Our framework is inspired by the kth order total generalized variation (TGV k ) [6,7], which is known to be effective in restoring the piecewise polynomial image with sharp edges [36]. Specifically, following the SLRM framework for the piecewise constant image [52], we assume that the image singularities (including both jumps and hidden jumps) are located in the zero level set of a band-limited periodic function (called the annihilating polynomial ). Then we can derive that the gradient can be decomposed into another vector field and the residual, and the Fourier samples of the residual and the symmetric gradient of this vector field can be annihilated by the convolution with the Fourier coefficients of the annihilating polynomial (called the annihilating filter ). From these annihilation relations, we deduce that the multi-fold Hankel matrices generated from the Fourier samples of a piecewise smooth image are low rank, which in turn enables a balance between the derivatives of order 1, . . . , k via a combined rank minimization of multi-fold Hankel matrices.
As a by-product of the proposed structured low rank matrix framework, we further introduce a wavelet frame based sparse regularization model for the piecewise image restoration via the continuous domain regularization. Briefly speaking, if we can associate a signal/image with a low rank Hankel matrix, its right singular vectors form tight frame filter banks under which the canonical coefficients have a group sparsity according to the index of filters [10]. Then motivated by [35], we assume that the right singular vectors are estimated from a pre-restoration process, and we propose a sparse regularization via the wavelet frame analysis approach (e.g. [13]) as an image restoration model via the continuous domain regularization. Notice that there are several wavelet frame based approaches for the piecewise smooth image restoration [12,24,41] in the literature. However, while these existing approaches focus on the sparse approximation of a discrete image, our approach comes from the relaxation of the structured low rank matrices generated by the Fourier samples for the continuous domain regularization.

SLRM for piecewise smooth image restoration
Our SLRM framework for the piecewise smooth function is mostly related to the recent extension of the SLRM framework to the two dimensional functions in [49,50,51,52,54]. Briefly speaking, if the singularity curves of a piecewise constant/holomorphic function, i.e. the supports of the first order (real/complex) derivatives of a target image, lie in the zero level set of an annihilating polynomial, the Fourier transform of the derivatives can be annihilated by the convolution with the annihilating filter. This annihilation relation in turn corresponds the Fourier samples of derivatives to the structured low rank matrices. Based on this framework, the SLRM framework for the piecewise constant image restoration is proposed and studied in [50,51,52,53], together with a restoration guarantee [49].
This annihilation relation of the gradient can be easily extended to the higher order derivatives whenever the jump discontinuities of the corresponding derivatives are located in the zero level set of a trigonometric polynomial [50]. For instance, we can derive an annihilation relation for the piecewise linear function by considering the Fourier transform of the second order derivatives. Based on this idea, the authors in [39] proposed a so-called generalized structured low rank (GSLR) approach for the piecewise smooth image restoration. Inspired by the 2-fold inf-convolution [20], the GSLR approach restores a superposition of a piecewise constant layer and a piecewise linear layer whose first and second order derivatives correspond to low rank Hankel matrices in the frequency domain respectively. By balancing the first order and second order derivatives, the GSLR approach has demonstrated significant improvements in the piecewise smooth image restoration tasks over the existing approaches.
Note that it is not difficult to extend the previous GSLR framework to generic piecewise smooth functions. More precisely, by considering the annihilation relation of higher order derivatives, we can extend the GSLR to the so-called k-fold inf-convolution (k ≥ 3) for functions with higher regularity. However, in a piecewise smooth function, there could exist singularities on which the derivatives have jump discontinuities. In this case, the annihilating polynomial describing the image singularities cannot annihilate the kth order derivatives [50,Proposition 4]. Since we then need to consider the power of the annihilating polynomial, or equivalently, increase the size of the filter to guarantee the annihilation relation for higher order derivatives, the GSLR framework can degrade the low rank (multi-fold) Hankel matrix structure of the higher order derivatives.
Unlike the GSLR framework in [39] based on the inf-convolution, the proposed SLRM framework is based on the total generalized variation framework [7]. Briefly speaking, we decompose the gradient into the residual whose Fourier samples correspond to the structured low rank matrix and a vector field whose Fourier samples of the symmetric gradient correspond to the structured low rank matrix. In other words, since our proposed SLRM framework uses the Fourier samples of the successive first order derivatives, it does not require the powers of the annihilating polynomial describing the image singularities, which enables us to obtain better low rank Hankel matrix structures corresponding to the piecewise smooth functions. In addition, since the TGV takes the inf-convolution as a special case, it can also be verified that the proposed SLRM framework is more general than the GSLR framework for the piecewise smooth functions.

Organization and notation of paper
The rest of this paper is organized as follows. In section 2, we present the proposed structured low rank matrix framework for the piecewise smooth functions. We first describe the one dimensional case to see the idea clearly, then we extend to the idea to the two dimensional framework. In section 3, we present an image restoration model based on the wavelet frame as an application of the proposed structured low rank matrix framework, followed by an alternating minimization algorithm, and some numerical results are presented to demonstrate the performance of our framework in the piecewise smooth image restoration in section 4. Finally, section 5 concludes this paper with a few future directions. All technical proofs will be postponed to appendices.
Throughout this paper, all two dimensional images and two dimensional k tensors defined on the discrete grid will be denoted by the bold faced lower case letters. Note that a two dimensional discrete image and a two dimensional discrete k tensor can also be identified with a vector and a 2 k -tuple of vectors (as well as a sequence or a 2 k -tuple of sequences supported on the grid) respectively whenever convenient. All matrices will be denoted by the bold faced upper case letters, and the mth row and the nth column of a matrix Z will be denoted by Z (m,:) and Z (:,n) , respectively. Denote by with N ∈ N, the set of N × N grid. The space of complex valued functions on O and the space of complex k-tensor valued functions on O are denoted by V C |O| and V k C |O|×2 k , respectively. Notice that V = V 0 . Given two rectangular grids K and M, we define Operators on both images and tensors are denoted as the bold faced caligraphic letters. For instance, letting v ∈ V and K be a rectangular K 1 × K 2 grid, the corresponding Hankel matrix Hv is an M 1 × M 2 matrix (M 1 = |O : K| and M 2 = |K|) generated by concatenating K 1 × K 2 patches of v into row vectors. Note that, in the sense of multi-indices, we have With a slight abuse of notation, we also use to denote the k-fold Hankel matrix constructed from q = (q 1 , . . . , q 2 k ) ∈ V k .

Structured low rank matrix framework for piecewise smooth functions
In this section, we introduce our structured low rank matrix framework for piecewise smooth functions. For simplicity, we consider the piecewise linear function throughout this paper. Note, however, it is not difficult to extend the proposed framework into the general piecewise smooth functions.

SLRM framework for one dimensional signals
We first establish the structured low rank matrix framework of the following one dimensional piecewise linear model from its Fourier sample In (2.1), α j , β j ∈ C, 1 [xj ,xj+1) denotes the characteristic function on the interval [x j , x j+1 ): , and 0 otherwise, and −1/2 < x 1 < x 2 < · · · < x K < 1/2 are the location of singularities, i.e. by the singularities we mean they include both the jumps and the hidden jumps (jumps of derivatives).
In the sense of distribution, the derivative u satisfies Hence, the Fourier transform of u − p is a linear combination of complex sinusoids: In addition, since the derivative of p is also a Dirac stream: its Fourier transform is expressed as i.e. a linear combination of complex sinusoids. We introduce the following trigonometric polynomial Since ϕ(x j ) = 0 for j = 1, . . . , K, it follows that ϕ (u − p) = 0 and ϕp = 0 in the sense of distribution. In the frequency domain, since we have it follows that for ξ ∈ R. Therefore, the Fourier transforms of both u − p and p are annihilated by the convolution under the Fourier coefficients of ϕ. In many practical cases, we consider the contiguous Fourier samples (i.e. the samples on a regular grid In the matrix-vector multiplication form, (2.10) and (2.11) lead to which shows that two (N −K+1)×K Hankel matrices H F (u − p) [N ] and H F (p ) [N ] have nontrivial nullspaces. In particular, for M > K, we can see that which shows that it is possible to convert the piecewise regularity in the continuous domain into the low rank Hankel matrices corresponding to the discrete Fourier samples.

SLRM framework for two dimensional images
Now we aim to establish the structured low rank matrix framework for the following two dimensional piecewise linear function model from its Fourier samples where the sample grid O is defined as (1.2). In (2.12), α j = (α j1 , α j2 ) ∈ C 2 , β j ∈ C and 1 Ωj denotes the characteristic function on a domain Ω j ; 1 Ωj (x) = 1 if x ∈ Ω j , and 0 otherwise. Without loss of generality, we assume that Ω j lies in [−1/2, 1/2) 2 for simplicity. (However, it is not difficult to generalize the setting into an arbitrary rectangular We further assume that (2.12) is expressed with the smallest number of characteristic functions such that Ω j 's are pairwise disjoint. Then the singularities of u, including both jump discontinuities of u and hidden jump discontinuities of u (jumps of derivatives), agree with Γ = J j=1 ∂Ω j , which will be called the singularity set of u throughout this paper. Generally, it is difficult to directly establish the SLRM framework without any further information on the singularity set Γ. Inspired by the two dimensional FRI framework for the piecewise constant function [50,51,52,53], we assume that there exists a finite rectangular and symmetric 4 index set K such that (2.14) Throughout this paper, we call any function ϕ(x) in the form of (2.14) the trigonometric polynomial, and the zero level set x ∈ R 2 : ϕ(x) = 0 the trigonometric curve. For a trigonometric polynomial ϕ in (2.14), the degree of ϕ is defined as an ordered pair of degrees in each coordinate, and is denoted by deg(ϕ). In particular, ϕ in (2.14) with the smallest degree is called the minimal polynomial. The algebraic properties of trigonometric polynomials and curves, including the existence of the minimal polynomial, are studied in [52], to which interested readers can refer for details. Under this setting, we present theorem 2.1 to establish the following (linear) annihilation relation for the piecewise linear function. The proof can be found in Appendix B.
Theorem 2.1. Let u(x) be defined as in (2.12) where the singularity set Γ satisfies (2.14). Then for some vector field p = (p 1 , p 2 ), we have and the Fourier transform F is performed to each component of ∇u − p and ∇ s p.
Based on theorem 2.1, we call ϕ satisfying (2.15) and (2.16) an annihilating polynomial (for u under a vector field p), and the Fourier coefficients a := {a(k) : k ∈ K} an annihilating filter. Note that the gradient of u is decomposed into ∇u − p and p such that the Fourier transforms of both ∇u − p and ∇ s p are annihilated by the same annihilating filter a. In addition, since the proof of theorem 2.1 tells us that p is the piecewise constant vector field in ∇u, (2.16) can be viewed as an extension of annihilation relation for piecewise constant functions to piecewise constant vector fields.
Remark 2.2. According to the proof, theorem 2.1 can also be written as follows: for a function u defined as in (2.12), ∇u can be decomposed into with a piecewise constant vector field p and a Radon vector measure ν supported on Γ such that We mention that the decomposition (2.18) is well-defined in the sense that F(∇ s ν) and F(∇ s p) do not share the minimal annihilating polynomial. To see this, let ϕ in (2.14) be the minimal polynomial for Γ. For Since ϕ is the minimal polynomial for a trigonometric curve Γ, ∇ϕ = 0 a.e. on Γ [52,Proposition A.4]. Hence, ϕ∇ s ν = 0, or equivalently, F(∇ s ν) does not satisfy the annihilation relation with a. In fact, ∇ s ν is annihilated by ϕ 2 , and since ϕ is the minimal polynomial for Γ, ϕ 2 is the minimal polynomial for annihilating ∇ s ν.
We claim that the linear annihilation relations (2.15) and (2.16) can be seen as balancing the annihilations of first and second order derivatives. To see this, we firstly consider an extreme case; a piecewise constant function case. In this case, it would be favorable to choose p = 0, so that (2.15) reduces to which is the annihilation relation for ∇u in [50,51,52,53]. For the second extreme case, we note that In addition, when ∇u is a piecewise constant vector field (i.e. u is continuous on Γ), we can simply choose which in turn measures the jump discontinuities of u and p in our setting. In particular, if u is continuous on ∂Ω j for some j, this ∂Ω j will not be reflected in (2.15). We also note that if we restrict p in the range of ∇, i.e. p = ∇u 2 , then by letting u 1 = u − u 2 , we can rewrite (2.15) and (2.16) as the GSLR framework in [39], which promotes a similar balance through a decomposition u = u 1 + u 2 with a piecewise constant u 1 and a piecewise linear u 2 . We can easily see that (2.24) and (2.25) are closely related to the following inf-convolution Noting that when p = ∇u 2 , then by letting u 1 = u − u 2 , the TGV (2.23) becomes the above inf-convolution, which means that the TGV takes the inf-convolution as a special case. Likewise, we can see that (2.15) and (2.16) in theorem 2.1 take (2.24) and (2.25) in [39] as a special case.
Remark 2.3. We further mention that theorem 2.1 is different from the GSLR framework in [39]. To see this, we assume that the singularity set Γ satisfies (2.14) with the minimal polynomial ϕ, i.e. K is the smallest support of the Fourier coefficients a. Since ∇ϕ = 0 a.e. on Γ [52, Proposition A.4], ϕ may not be able to annihilate the second order derivatives of a piecewise linear function in general [50]. Hence, (2.24) and (2.25) prefer to decompose u = u 1 + u 2 where u 1 is a piecewise constant function and u 2 is a piecewise linear spline (a continuous piecewise linear function). Meanwhile, the proposed framework is established by annihilating the jumps for each derivative successively, which enables to cover a broader range of piecewise linear functions. Hence, our framework is more generalized than the GSLR in [39].
Next, we will present two examples to illustrate this.
Then we can see that is a minimal polynomial for the singularity set {−1/4, 0, 1/4}. Since Hence, it follows that In other words, both frameworks can establish annihilation relations if a given function can be decomposed into a piecewise constant function and a piecewise linear spline.
being a minimal polynomial for the singularity set {−1/4, 1/4}. In this case, we have If we choose p = 1 [−1/4,1/4) , then we have However, we cannot find a piecewise constant u 1 and a piecewise linear spline u 2 such that Notice that we have Hence, these two examples again illustrate that the proposed framework is more generalized than the GSLR framework in [39].
In our setting, the Fourier transform of u is sampled on the grid O in (1.2) with N ∈ N large enough to guarantee a high image resolution. Hence, (2.15) and (2.16) become the following finite systems of linear equations where m ∈ O : K. In the matrix-vector multiplication form, we have O have nontrivial nullspaces. In addition, when the filter support K defining two multi-fold Hankel matrices is sufficiently large, both of them have nontrivial nullspaces as well. To see this, let ϕ be the minimal polynomial for the singularity set with coefficients a supported on K. Then for any trigonometric polynomial η with coefficients b supported on K : K such that a * b is supported on K , it follows that is an annihilating polynomial and a * b is an annihilating filter, due to the associativity of convolution. In particular, by letting we can see that is also an annihilating polynomial, or equivalently, the translation a(· − m) also satisfies (2.29) and (2.30).
Based on this observation, we present proposition 2.6 to demonstrate the low rank properties of multi-fold Hankel matrices, which will establish the correspondence between the Hankel matrices and the complexity of the singularity set of u.
Proposition 2.6. Let u(x) be defined as in (2.12) with the singularity set Γ = J j=1 ∂Ω j satisfying (2.14), where ϕ is the minimal polynomial with coefficients on K. For an assumed filter support K strictly containing K, we have O are of low-rank, which enables to convert the piecewise regularity in the continuous domain into the low rank multi-fold Hankel matrices corresponding to the discrete Fourier samples. Notice that the GSLR framework in [39] promotes the decomposition u = u 1 + u 2 where the Fourier samples of ∇u 1 and ∇ 2 u 2 correspond to the low rank Hankel matrices. In contrast, the proposed SLRM framework is established by decomposing the gradient of u into ∇u − p and p such that the Fourier samples of ∇u − p and ∇ s p correspond to the low rank Hankel matrices. See figure 1 for the schematic illustrations.

Continuous domain regularization for piecewise smooth image restoration
Let f ∈ V be a degraded measurement modeled as

1)
Fourier sample Low rank Hankel matrix Low rank Hankel matrix  where v = F(u) O with u defined as in (2.12), and ζ is some measurement error 5 . According to section 2, the (multi-fold) Hankel matrices corresponding to are low rank. Hence, we can consider to restore v ∈ V, as a continuous domain regularization for the piecewise smooth image restoration. Here, for k = (k 1 , k 2 ) ∈ O, respectively. When q = 0, (3.2) reduces to which is a rank minimization model based on the GSLR framework. From the two extreme cases (3.5) and (3.6), we can see that (3.2) aims to restore piecewise smooth functions by decomposing Dv into Dv − q and q, thereby balancing the low rank multi-fold Hankel matrices of the Fourier samples of the first order derivatives and the second order derivatives. In addition, we can also see that (3.2) takes the GSLR rank minimization model (3.8), which directly decomposes v = v 1 + v 2 , as a special case. Since our SLRM framework is more generalized than the GSLR framework (see remark 2.3 and examples 2.4 and 2.5), it can be expected that (3.2) is able to restore a wider range of piecewise smooth functions than (3.8).

From low rank model to tight frame approach
Though (3.2) is an NP-hard problem, there are numerous tractable approaches available, including the convex nuclear norm relaxation (e.g. [14,30]), the iterative reweighted least squares (IRLS) for the Schatten p-norm minimization [33,39,48,53], etc. Inspired by the SVD of a low rank Hankel matrix, we propose another relaxation of (3.2), similar to [10,Theorem 3.2] where the continuous domain regularization for the piecewise constant image restoration is studied. For this purpose, we present the main idea of the image restoration model in theorem 3.1. The proof is postponed to Appendix C.
We begin with introducing some notation. For a set of K 1 × K 2 filters a 1 , . . . , a M2 supported on K, we define W and W * (the adjoint of W) as where S a is a discrete convolution under the periodic boundary condition: In other words, both W and W * are concatenations of discrete convolutions. by reformulating each column vector into a K 1 × K 2 filter supported on K. Then W 1 and W 2 defined as (3.9) by using filters {a 11 , . . . , a 1M2 } and {a 21 , . . . , a 2M2 } satisfies

15)
where the discrete convolution is performed on each component of Dv−q and Eq, respectively. Consequently, if H (Dv − q) (and H (Eq), respectively) is of low rank, then its right singular vectors construct a tight frame under which Dv − q (and Eq, respectively) is sparsely represented.
In words, theorem 3.1 tells us that if the SVDs of multi-fold Hankel matrices are known as oracles, we can explicitly construct tight frames under which Dv − q and Eq are sparsely represented, and the sparsity of canonical coefficients can be grouped according to the filters. Hence, motivated by the idea in [35], assume that v ∈ V and q ∈ V 1 are the a-priori estimations of v and q with the SVDs Then we define the tight frame transforms W 1 and W 2 in (3.9) via a 1l = M Since (3.14) and (3.15) are then approximately true under these W 1 and W 2 , we remove the group sparsity pattern in the canonical coefficients for the better sparse approximation instead. This leads us to solve where the 1 norms take the form of with some ν k > 0 and a small ε > 0 to avoid the division by zero. Hence, we relax the sparsity of tight frame coefficients over the entire range of W k 's (not necessarily in groups as in (3.14) and (3.15)), we expect to achieve more flexibility, thereby leading to the improvements in restoration performance.
Since the wavelet frame based relaxation model (3.16) is inspired by the (3.2) via the SVDs of H (Dv − q) and H (Eq), we can say that our relaxation model (3.16) restores piecewise smooth functions by decomposing Dv into Dv − q and q, thereby balancing the low rank Hankel matrices of the Fourier samples of the first order derivatives and that of the second order derivatives. In addition, we further note that, following [10], it is also possible to consider the following data driven tight frame model with the 0 norm γ k · c k 0 (k = 1, 2) encoding the number of nonzero entries in c k 's, to learn tight frames and restore the Fourier sample v simultaneously. Even though it is not clear at this point whether the adaptive tight frame system will lead to better restoration results or not, throughout this paper, we only consider the model (3.16) rather than the data driven tight frame model (3.17), and the reasons are as follows. First of all, given that W k 's are properly estimated, it may not be necessary to further learn them with additional computational costs. Second, since (3.16) is convex whereas (3.17) is nonconvex, we can easily expect better behavior and theoretical support for the numerical algorithms. Most importantly, from the viewpoint of wavelet frame based image restoration, the model (3.16) is an analysis approach [13] while the data driven tight frame model (3.17) can be classified into the balanced approach [9,21]. (See Appendix A for the brief descriptions on the wavelet frame based models.) It is well known that the analysis approach reflects the structure of a target image better than other approaches (e.g. [11]). Since we derive the wavelet frame based approach from the structured low rank matrix framework, it can be expected that the convex relaxation model (3.16) will reflect our SLRM frameworks for the piecewise smooth functions better than the data driven tight frame model (3.17).

Alternating minimization algorithm
Among numerous algorithms which can solve the convex model (3.16), we adopt the ADMM [28] or the split Bregman algorithm [36], which can convert (3.16) into several subproblems with closed form solutions, together with the convergence guarantee [13]. More precisely, let c 1 = W 1 (Dv − q), and c 2 = W 2 (Eq).
Under this reformulation, the overall algorithm is summarized in algorithm 1. For (3.18), since W * k W k = I for k = 1, 2, we solve the following system of linear equations Depending on the formulation of A, various methods can be used to solve (3.23) efficiently. For example, when A is a pointwise multiplication in the frequency domain (e.g. image denoising, image deblurring, and CS restoration), so are the constituent operators K 11 , K 21 , and K 22 , and thus we can solve (3.23) by the pointwise Cramer's rule. For a more generalized A, we can use a distributed optimization based method [5]. The closed form solutions for (3.19) and (3.20) is expressed in terms of the soft thresholding: More precisely, the soft thresholding operator T γ (c) for c ∈ V M2 k and γ = γ 1 · · · γ M2 T is defined as the following componentwise manner: for k ∈ O, l = 1, . . . , M 2 , and m = 1, . . . , 2 k , with the convention that 0/0 = 0.

Numerical results
In this section, we conduct some numerical simulations in the context of restoration from the partial Fourier samples as a proof-of-concept study. Specifically, to compare the performance of the direct rank minimization (3.2) and the wavelet frame relaxation (3.16), we choose to compare the relaxation model solved by algorithm 1 with the following Schatten 0-norm relaxation [33,39,48,53] (SLRM model) of (3.2): where R M denotes a restriction onto the known sample grid M. In addition, since it has been demonstrated in [39] that the GSLR framework outperforms the structured low rank matrix approaches which consider either the first order derivatives or the second derivatives, we also compare with the following GSLR model in [39]: In (4.2) and (4.3), Z 0 is the Schatten 0-norm of a matrix Z defined as Z 0 = ln det (Z * Z) 1/2 + εI with a small constant ε > 0, and (4.2) and (4.3) are solved by the iterative reweighted least squares method described in [39]. Finally, to further study the improvements over the conventional on-the-grid approaches, we compare with the piecewise linear framelet (Fra) model (e.g. [13]) the TGV model [7,44] min u,p 5) and the inf-convolution (IF) model [20] min u1,u2 where F denotes the two dimensional discrete Fourier transform. Throughout this paper, we use the split Bregman algorithm (e.g. [13,34,36]) to solve (4.4)-(4.6). All experiments are implemented on MATLAB R2014a running on a laptop with 64GB RAM and Intel(R) Core(TM) CPU i7-8750H at 2.20GHz with 6 cores. Throughout the experiments, we test two synthetic images ("Ellipses" and "Rectangles"), and two natural images ("Airplane", and "Car") taking the values in [0, 1], as shown in figure 2. The data f is synthesized by randomly sampling 20% of 256 × 256 Fourier samples via the variable density sampling method described in [47]. A complex white Gaussian noise with the standard deviation 1 is also added to generate a noisy partial sampling. For (4.1)-(4.3), we use the K × K square patch for simplicity. Specifically, we choose K = 31 for the "Ellipses", K = 25 for the "Rectangles", and K = 51 for "Airplane" and "Car", depending on the geometry of the target images. More precisely, in (4.1)-(4.3), we use the discrete Fourier transform of the images restored by the TGV model (4.5) and the inf-convolution model (4.6) respectively, to compute the SVDs of multi-fold Hankel matrices using the K × K patches, which will be used for W k 's in (4.1) and initializations in (4.2) and (4.3). For (4.4), we choose W to be the undecimated tensor product piecewise linear B-spline framelet transform with 1 level of decomposition [26]. For (4.5) and (4.6), we use the forward difference with the periodic boundary condition for the difference operators ∇, ∇ s , and ∇ 2 . In all of the experiments, we have manually tuned the regularization parameters to achieve the optimal restoration results in each scenario. For the quantitative comparison, we compute the signal-to-noise ratio (SNR), the high frequency error norm (HFEN) [59], and the structure similarity index map (SSIM) [64]. Note that for (4.1)-(4.3), the restored image is computed via the inverse DFT of the restored Fourier samples. table 1 summarizes the SNR, the HFEN, and the SSIM of the aforementioned approaches. figures 3 to 6 display the visual comparisons (the first row) together with the zoom-in views (the second row) and the error maps (the third row) of (4.1) against (4.2)-(4.6). Throughout this paper, all restored images are displayed in the window level [0, 1], and all error maps are displayed in the window level [0, 0.2] for fair comparisons. We can see that the relaxed model (4.1) consistently performs best, and the SLRM model (4.2) performs the second in almost every index in all scenarios, and the improvements are visually observable as well, both of which demonstrate that our proposed SLRM framework performs well in the piecewise smooth image restoration.
At first glance, since the models (4.1) and (4.2) can be regarded as different relaxations of the structured low rank matrix framework for the continuous domain regularization, the results also suggest that the off-the-grid regularization performs better than the on-the-grid regularization as it can reduce the basis mismatch between the true support (or the true singularity) in continuum and the discrete grid, leading to the improvements in the indices. In fact, due to such a basis mismatch, the conventional on-the-grid approaches ((4.4)-(4.6)) suffers from the errors concentrated near the image singularities as well as the distorted shapes compared to (4.1), which can be seen in the error maps and the zoom-in views of figures 3 to 6. Most importantly, the numerical results illustrate the followings. First of all, from the comparisons between (4.2) and (4.3), we can see that the proposed SLRM framework performs better than the GSLR framework in [39] for the piecewise smooth image restoration. Second, from the comparisons between (4.1) and (4.2), we can further see that the relaxation into the wavelet frame analysis approach can achieve further improvements over the direct rank minimization. As previously mentioned, since the proposed SLRM framework takes the GSLR framework in [39] as a special case, we can restore broader range of piecewise smooth images. In addition, it should be noted that the annihilation of the Fourier samples of the second order derivatives in general require larger annihilating filter than the first order derivatives [50], which may degrade the low rank structures of the corresponding multi-fold Hankel matrix. As a consequence, the GSLR framework can introduce staircase artifacts in the smooth region compared to the proposed SLRM framework, as can be seen in the zoom-in views of figures 3 to 6. Finally, since the weights corresponding to the derivatives in the image domain amplifies the noise in the high frequencies, it is likely that the direct rank minimization approaches ((4.2) and (4.3)) lead to the artifacts near the image singularities corresponding to the high frequency components, whereas the relaxation into the wavelet frame analysis approach ((4.1)) can achieve the denoising effect in spite of the amplified noise, leading to the better restoration results with less artifacts near the edges, as can be seen in the error maps in figures 3 to 6.

Conclusion and future directions
In this paper, we have introduced a new structured low rank matrix framework for the piecewise smooth image restoration. Following the previous structured low rank matrix framework for the piecewise constant image [52], we assume that the image singularities lie in the zero level set of a band-limited periodic function. Inspired by the total generalized variation [7], we derive the annihilation relation of the Fourier samples of the gradient and the (symmetric) gradient tensor, which in turn leads to the low rank multi-fold Hankel matrices to balance between the first order derivatives and the second order derivatives. In addition, as a by-product of the proposed structured low rank matrix framework, we further introduce a wavelet frame based sparse regularization model for the piecewise image restoration via the continuous domain regularization, from the SVDs of low rank multi-fold Hankel matrices. Finally, the numerical experiments show that the proposed wavelet frame based model based on the proposed SLRM framework outperforms both the conventional on-the-grid approaches and the existing the SLRM framework as well as the existing rank minimization approaches for the piecewise smooth image restoration. For the future work, we plan to generalize the frameworks in [52] to the piecewise smooth image restoration. More precisely, we need to develop sampling guarantees for the unique restoration of a continuous domain piecewise smooth image from its uniform low pass Fourier samples, based on the proposed annihilation relations. Unfortunately, it is not clear for us at this moment under what conditions we can restore the singularity set from uniform Fourier samples as the annihilation relations involve a blind deconvolution problem due to the additional unknown vector field p in theorem 2.1. Nevertheless, this is definitely a future direction we would like to work on. We are also interested in proving the restoration guarantee of (3.2) (or its convex nuclear norm relaxation) given that A = R M , i.e. the sample grid M is drawn from the image grid O uniformly at random, which generalizes the theoretical frameworks in [49] to the piecewise smooth image restoration.
Finally, to broaden the scope of applications, it is also likely to apply the idea in this paper to the various image restoration tasks such as the blind deconvolution (especially that in the frequency domain which is related to the retinex [46]) and the sparse angle CT. In fact, it is not clear at this moment how to extend this approach to the sparse angle CT restoration, as the measurement can be viewed as the Fourier samples on the discrete radial grid. Nevertheless, this is also an important future work for the broader scope of applications. Finally, observing the successive structure in our SLRM framework, we are also interested in developing a deep learning framework, motivated by the recent works on the deep learning techniques combined with the structured low rank matrix approaches [40,43,56,57].

Appendix A. Preliminaries on tight wavelet frames
Provided here is a brief introduction on tight wavelet frames and the wavelet frame based image restoration. Interested readers may consult [9,26,27,60,62] for detailed surveys on tight wavelet frames, and [11] for detailed reviews on the wavelet frame based image restoration. For the sake of simplicity, we only discuss the real valued wavelet tight frame systems, but note that it is not difficult to extend the idea to the complex case. Let H be a Hilbert space equipped with an inner product ·, · . A (Bessel) sequence {ϕ n : n ∈ Z} ⊆ H is called a tight frame on H if Given {ϕ n : n ∈ Z} ⊆ H, we define the analysis operator W : H → 2 (Z) as u ∈ H → Wu = { u, ϕ n : n ∈ Z} ∈ 2 (Z).
The synthesis operator W T : 2 (Z) → H is defined as the adjoint of W: Then {ϕ n : n ∈ Z} is a tight frame on H if and only if W T W = I. It follows that, for a given tight frame {ϕ n : n ∈ Z}, we have the following canonical expression: with Wu = { u, ϕ n : n ∈ Z} being called the canonical tight frame coefficients. Hence, the tight frames are extensions of orthonormal bases to the redundant systems. In fact, a tight frame is an orthonormal basis if and only if ϕ n = 1 for all n ∈ Z.
One of the most widely used class of tight frames is the discrete wavelet frame generated by a set of finitely supported filters {a 1 , · · · , a m }. Throughout this paper, we only discuss the two dimensional undecimated wavelet frames on 2 (Z 2 ), but note that it is not difficult to extend to 2 (Z d ) with d ≥ 3. For a ∈ 1 (Z 2 ), define a convolution operator S a : 2 (Z 2 ) → 2 (Z 2 ) by Given a set of finitely supported filters {a 1 , · · · , a m }, define the analysis operator W and the synthesis operator W T by respectively. Then, the direct computation can show that the rows of W form a tight frame on 2 (Z 2 ) (i.e. W T W = I) if and only if the filters {a 1 , · · · , a m } satisfy m j=1 l∈Z 2 called the unitary extension principle (UEP) condition [37]. Finally, for a two dimensional discrete image on the finite grid, S a in (A.2) with a finitely supported a denotes the discrete convolution under the periodic boundary condition throughout this paper. For the image restoration, we assume that there exists a tight wavelet frame W defined as (A.3) with filters {a 1 , · · · , a m } such that u is sparse under W. This leads us to solve and restore u = W T c. In (A.6), the first term is a data fidelity term, the third 1 norm is defined as with the parameters γ = [γ 1 , . . . , γ m ] T and γ j > 0, to promote the sparsity of c. The second term penalizes the distance between c and its projection onto R(W). Since we obtain the image via u = W T c, the second term in fact forces the sparse coefficient c close to the canonical coefficient of the restored image. Since the magnitude/decay of the canonical coefficient reflects the regularity of the image under some mild conditions on the tight frame system W [4], we can see that the second and third terms promote a balance between the sparsity of coefficient and the regularity of image. We can further impose the flexibility on (A.6) by introducing the weight on the second term: called the balanced approach [9,21]. If µ = 0, then (A.8) becomes the synthesis approach [25,29,31] min as we aim to find the sparsest coefficient that synthesizes the image. On the other hand, if µ = ∞, then it must be (I − WW T )c = 0. Since we then have c ∈ R(W), it follows that c = Wu for some u. Hence, (A.6) becomes the analysis approach [13] min u 1 2 Since the analysis approach promotes the sparse canonical coefficient, it emphasizes the regularity of the restored image. These three approaches are equivalent if and only if W T W = WW T = I, i.e. W forms an orthonormal basis. However, since we consider the redundant system, each approach will lead to the different restoration results. Notice that the coefficient vector c obtained from the synthesis approach (A.9) is much sparser than the canonical coefficient Wu obtained from the analysis approach (A.10). It is empirically observed that the synthesis approach (A.9) tends to yield some artifacts in the restored image. In contrast, the analysis approach (A.10) usually have less artifacts as it promotes the regularity of the image along the singularities. Theoretically, it has been proved in [11] that, under a suitable choice of parameters, the analysis approach (A.10) can be seen as sophisticated discretization of a variational model including the Rudin-Osher-Fatemi model [61], which enables a geometric interpretation on (A.10).
Then it is obvious that we have Since |n j | = 1 almost everywhere on ∂Ω j , α T j x + β j is continuous on a compact set ∂Ω j , ∇u − p in (B.2) defines a finite Radon vector measure on R 2 supported on Γ. In addition, since (B.1) as well as (B.2) holds in the sense of tempered distribution, we can compute F(∇u − p)(ξ) as a Fourier transform of a measure (e.g. [32]). Namely, we have Let ϕ ∈ S, where S is the space of Schwartz functions (e.g. [32]). Then the direct computation shows that (F(∇u − p) * ϕ) (ξ) satisfies (F(∇u − p) * ϕ) (ξ) = This proves (2.15).
To complete the proof, we further note that in the sense of distribution, ∇ s p satisfies Following the similar arguments, it is not difficult to verify that dΦ j in (B.5) defines a finite Radon tensor measure on R 2 supported on Γ. which means that F(∇ s p) is defined as This means that W 1 and W * 1 defined as (3.9) and (3.10) using a 11 , . . . , a 1M2 satisfies S a 1l S a 1l (−·) (Dv − q) = Dv − q, which shows that the filters a 11 , . . . , a 1M2 form a tight frame system. In addition, we note that ∈ C M2 into a K 1 × K 2 filter supported on K. Then (C.1) implies that for l = r 1 + 1, . . . , M 2 , we have S a 1l (−·) (Dv − q) (k) = 0, k ∈ O : K, which proves (3.12) and (3.14). Similarly, we can prove that a 21 , . . . , a 2M2 satisfies M2 l=1 m∈K a 2l (k + m)a 2l (m) = δ(k), so that S a 2l S a 2l (−·) (Eq) = Eq.
This completes the proof.