---
abstract: 'Recently, the finite-rate-of-innovation (FRI) based continuous domain restoration is emerging as an alternative to the conventional on-the-grid compressed sensing (CS) MRI restoration due to its ability to alleviate the basis mismatch between the true support of the shape in the continuous domain and the discrete grid. In this paper, we propose a new off-the-grid approach for the CS-MRI reconstruction. Following the recent works on two dimensional FRI, we assume that the discontinuities/edges of the image are localized in the zero level set of a band-limited periodic function. This assumption induces the linear dependencies among the Fourier samples of the gradient of the image, which leads to a low rank two-fold Hankel matrix. We further observe that the singular value decomposition of a low rank Hankel matrix corresponds to an adaptive tight frame system which can represent the image with a few nonzero canonical coefficients. Based on this observation, we propose a data driven tight frame based off-the-grid CS-MRI reconstruction model for the MR image reconstruction. To solve the nonconvex and nonsmooth model, a proximal alternating minimization algorithm with a guaranteed global convergence is adopted. Finally, the numerical experiments show that our proposed data driven tight frame based approach outperforms the existing approaches.'
author:
- 'Jian-Feng Cai[^1]'
- 'Jae Kyu Choi[^2]'
- 'Ke Wei[^3]'
title: '[[Off-the-Grid Compressed Sensing MRI Reconstruction via Data Driven Tight Frame Regularization]{}]{}[^4]'
---

Magnetic resonance imaging, finite-rate-of-innovation, structured low rank matrix completion, (tight) wavelet frames, data driven tight frames, proximal alternating schemes

42B05, 65K15, 65R32, 68U10, 90C90, 92C55, 94A12, 94A20

Introduction {#Introduction}
============

Magnetic resonance imaging (MRI) is one of the most widely used medical imaging modality in clinical diagnosis [@E.M.Haacke1999]. It is non-radioactive, non-invasive, and has excellent soft tissue contrasts such as T1 and T2 with high spatial resolution [@Y.Liu2015]. Among these merits, the availability of high spatial resolution images, which will be the focus of this paper, facilitates early diagnosis by enabling the detection and characterization of clinically important lesions [@G.Ongie2015; @E.VanReeth2012]. However, since the so-called k-space data acquisition is limited due to physical (gradient amplitude and slew-rate) and physiological (nerve stimulation) constraints [@C.M.Hyun2018; @Y.Liu2015; @M.Lustig2007], there has been increasing demand for methods which can reduce the amount of acquired data without degrading the image quality [@M.Lustig2007].

When the k-space data is undersampled, the Nyquist sampling criterion is violated, and this inevitably leads to the aliasing in the reconstructed image [@M.Lustig2007]. In the literature, the famous compressed sensing (CS) MRI can be viewed as a sub-Nyquist sampling method which exploits the sparse representation of an image to compensate the undersampled k-space data [@E.J.Candes2006; @Donoho2006; @C.M.Hyun2018; @M.Lustig2007]. A typical $\ell_1$ norm based CS-MRI model restores the MR image $\bu$ on the grid $\OO\subseteq\Z^2$ via $$\begin{aligned}
\label{CSType}
\min_{\bu}~\f{1}{2}\left\|\bmP_{\MM}\bmsF\bu-\bsf\right\|_2^2+\lambda\left\|\Ph\bu\right\|_1,\end{aligned}$$ where $\bmsF$ is the discrete Fourier transform (DFT), $\bsf$ is the undersampled k-space data on $\MM\subseteq\OO$, and $\bmP_{\MM}$ is the projection onto $\MM$. In addition, $\Ph$ is some sparsifying linear transformation, for example, by computing the discrete gradient [@M.Lustig2007] or the wavelet frame coefficients [@Y.Liu2015] of an image. The “on-the-grid” restoration model is widely used for regularizing smooth image components while preserving image singularities such as edges, ridges, and corners.

Even though the CS-MRI model \[CSType\] has shown strong ability to reduce the data acquisition time and thus received a lot of attention over the past few years [@Y.Liu2015], it has to be further improved as the k-space data is a discrete (and truncated) sampling of a Fourier transform of an underlying function in the *continuous domain* [@E.M.Haacke1999]. This means that, \[CSType\] will work well when the singularities of the image $\bu$ are well aligned with the grid $\OO$. However, even in the case of piecewise constant function whose (distribution) gradient is sparse in the continuous domain, its singularities may not necessarily agree with the discrete grid $\OO$, leading to the problem of *basis mismatch* [@Y.Chi2010; @G.Ongie2016; @J.Ying2017]. Such a basis mismatch between the true singularities and the discrete grid may result in the loss of sparsity of $\Ph\bu$, and thus degrades the performance of \[CSType\] [@J.Ying2017].

The direct restoration of continuous domain signal is emerging as a powerful alternative to the discrete domain CS approach [@B.N.Bhaskar2013; @E.J.Candes2014; @Y.Chen2014; @G.Ongie2018]. The ability of these “off-the-grid” schemes to allieviate the basis mismatch makes them especially attractive in signal restorations from partial Fourier samples [@G.Ongie2018]. One of the most successful examples is the finite-rate-of-innovation (FRI) framework [@T.Blu2008] which extends the Prony’s method [@Kay1988; @S.M.Kay1981; @Prony1795; @P.Stoica1997] to the so-called nonuniform “Dirac streams” with unknown knot locations. In addition, the recovery of Dirac streams from a few Fourier samples can be achieved by the atomic norm minimization [@E.J.Candes2013; @E.J.Candes2014; @G.Tang2015], which can be viewed as the continuous domain generalization of the CS framework. While these frameworks work well for one dimensional piecewise constant signals [@E.J.Candes2013] and two dimensional Dirac streams [@W.Xu2014], the extension to the (piecewise smooth) image restoration is not straightforward [@G.Ongie2016]. This is because the image singularities such as the edges and ridges in general form a continuous curve in a two dimensional domain, which violates the separation condition [@E.J.Candes2014] for the recovery guarantee.

Recently, there are several extensions of the FRI framework to the two dimensional continuous domain image restoration. Such extensions include the piecewise holomorphic complex image restoration [@H.Pan2014] and the piecewise constant real image restorations [@G.Ongie2018; @G.Ongie2015a; @G.Ongie2015; @G.Ongie2016]. These approaches are commonly based on the assumption that the singularity curves, i.e. the supports of the (real/complex) derivatives of a target image, lie in the zero level set of a band-limited periodic function, called the *annihilating polynomial*. The annihilation relation of the Fourier transform of the derivatives by convolution with the Fourier coefficients of the annihilating polynomial is then derived to recover the image singularities and eventually the image [@G.Ongie2016]. In addition, the authors in [@G.Ongie2016] established the necessary and sufficient conditions for the restoration of the edge set and the edge set aware recovery of the piecewise constant images from a finite number of uniformly sampled low frequency k-space data.

Given that the edges lie in the zero level set of an annihilating polynomial, the continuous domain image restoration is based on the fact that the structured matrix (Hankel/Toeplitz matrix) constructed from the Fourier samples of piecewise constant images is low rank [@G.Ongie2018; @G.Ongie2015; @G.Ongie2016]. Hence, we can apply the following structured low rank matrix completion model [@Y.Chen2014; @G.Ongie2017] $$\begin{aligned}
\label{LRSMCompletion}
\min~~\rank\left(\bmH\left(\La\bsv\right)\right)~~~~~~~\text{subject to}~~\bmP_{\MM}\bsv=\bsf\end{aligned}$$ to restore the fully sampled k-space data (the discrete sample of the Fourier transform of a piecewise constant function). Here, $\bmH$ is an operator that constructs a so-called two-fold Hankel matrix, and $\La$ is a weight matrix derived from the Fourier transform of $\na$ (the detailed definitions will be postponed until \[Preliminaries\]). By restoring the k-space data first, we can directly utilize the property of the underlying function in the continuous domain, thus achieving the better restoration results over \[CSType\] by reducing the basis mismatch.

Though \[LRSMCompletion\] is an NP-hard problem in general [@E.J.Candes2009], numerous tractable relaxation approaches have been proposed in the literature. One of them is the convex nuclear norm relaxation method (e.g. [@J.F.Cai2016; @M.Fazel2013]), together with the theoretical restoration guarantees (e.g. [@G.Ongie2018]). In addition, the iterative reweighted least squares (IRLS) approaches for the Schatten $p$-norm minimization are proposed in [@M.Fornasier2011; @K.Mohan2012; @G.Ongie2017], which can avoid the high computational cost of the SVD related to the rank minmization and the convex nuclear norm relaxation. Apart from these approaches, the nonconvex alternating projection methods based on the different parametrizations of the underlying low rank matrix structure are proposed and studied in [@J.F.Cai2018a; @J.F.Cai2019]. These nonconvex methods are reported to be superior to the other relaxation methods in terms of the computational efficiency while theoretical restoration guarantees are still available.

In this paper, we propose a new off-the-grid approach for the CS-MRI reconstruction based on the data driven tight frame [@J.F.Cai2014]. To this end, we first note that the SVD of a Hankel matrix induces a tight frame filter bank due to its underlying convolutional structure. More precisely, if we can associate a signal with a low rank Hankel matrix, its right singular vectors form a tight frame filter bank which allows us to represent the signal with a small number of nonzero canonical coefficients. Based on this observation, we develop the sparse regularization model via data driven tight frames [@J.F.Cai2014] for the CS-MRI restoration problem. Even though the data driven tight frame approach has been discussed for the sparse approximation of an image [@J.F.Cai2014; @J.K.Choi2018; @J.Wang2015; @R.Zhan2016], to the best of our knowledge, this is the first work that applies this approach for the MRI restoration based on a continuous image model. Finally, the numerical experiments show that our data driven tight frame approach outperforms the structured low rank matrix approaches based on \[LRSMCompletion\], leading to the state-of-the-art performance.

The rest of this paper is organized as follows. We first briefly review the related concepts, including (data driven) tight frames and structured low rank matrix approaches for CS-MRI, in \[Preliminaries\]. Then we present the off-the-grid CS-MRI reconstruction approach based on data driven tight frames and the alternating minimization algorithm in \[ModelandAlgorithm\], followed by the convergence of the algorithm. In \[Experiments\], experimental results are reported to demonstrate the performance of our new CS-MRI restoration method, and \[Conclusion\] concludes this paper with a few future directions.

Preliminaries and related works {#Preliminaries}
===============================

Introduction on tight frames {#TightFrames}
----------------------------

Here we provide a brief introduction on tight wavelet frames and data driven tight frames. Interested readers may consult [@J.F.Cai2008; @B.Dong2013; @B.Dong2015; @A.Ron1997; @Shen2010] for details on tight wavelet frames, and [@C.Bao2015; @J.F.Cai2014] for details on data driven tight frames. 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.

Denote by $\msH$ a Hilbert space and let $\la\cdot,\cdot\ra$ be the inner product defined on $\msH$. A countable set $\{\bvphi_n:n\in\Z\}\subseteq\msH$ is called a tight frame on $\msH$ if $$\begin{aligned}
\label{TightFrame}
\|\bu\|^2=\sum_{n\in\Z}|\la\bu,\bvphi_n\ra|^2~~~~~\text{for all}~~~\bu\in\msH.\end{aligned}$$ Given $\{\bvphi_n:n\in\Z\}\subseteq\msH$, we define the analysis operator $\bsW:\msH\to\ell_2(\Z)$ as $$\begin{aligned}
\bu\in\msH\mapsto \bsW\bu=\{\la\bu,\bvphi_n\ra:n\in\Z\}\in\ell_2(\Z).\end{aligned}$$ The synthesis operator $\bsW^T:\ell_2(\Z)\to\msH$ is defined as the adjoint of $\bsW$: $$\begin{aligned}
\bc\in\ell_2(\Z)\mapsto\bsW^T\bc=\sum_{n\in\Z}\bc[n]\bvphi_n\in\msH.\end{aligned}$$ Then $\{\bvphi_n:n\in\Z\}$ is a tight frame on $\msH$ if and only if $\bsW^T\bsW=\bsI$ where $\bsI$ is the identity on $\msH$. It follows that, for a given tight frame $\{\bvphi_n:n\in\Z\}$, we have the following canonical expression: $$\begin{aligned}
\bu=\sum_{n\in\Z}\la\bu,\bvphi_n\ra\bvphi_n,\end{aligned}$$ with $\bsW\bu=\{\la\bu,\bvphi_n\ra:n\in\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 $\|\bvphi_n\|=1$ for all $n\in\Z$.

One of the most widely used class of tight frames is the discrete wavelet frame generated by a set of finitely supported filters $\{\bq_1,\cdots,\bq_m\}$. In this paper, we only discuss the undecimated wavelet frames, which are also known as the translation invariant wavelet frame. For $\bq\in\ell_1(\Z)$, define a convolution operator $\bsS_{\bq}:\ell_2(\Z)\to\ell_2(\Z)$ by $$\begin{aligned}
\label{DiscreteConv}
(\bsS_{\bq}\bu)[n]=(\bq\ast\bu)[n]=\sum_{k\in\Z}\bq[n-k]\bu[k]~~~~~\text{for}~~~\bu\in\ell_2(\Z).\end{aligned}$$ Given a set of finitely supported filters $\{\bq_1,\cdots,\bq_m\}$, define the analysis operator $\bsW$ and the synthesis operator $\bsW^T$ by $$\begin{aligned}
\bsW&=\left[\bsS_{\bq_1[-\cdot]}^T,\bsS_{\bq_2[-\cdot]}^T,\cdots,\bsS_{\bq_m[-\cdot]}^T\right]^T,\label{AnalConv}\\
\bsW^T&=\left[\bsS_{\bq_1},\bsS_{\bq_2},\cdots,\bsS_{\bq_m}\right],\label{SyntConv}\end{aligned}$$ respectively. Then, the direct computation can show that the rows of $\bsW$ form a tight frame on $\ell_2(\Z)$ (i.e. $\bsW^T\bsW=\bsI$) if and only if the filters $\{\bq_1,\cdots,\bq_m\}$ satisfy the following [*unitary extension principle*]{} (UEP) condition [@B.Han2011]: $$\begin{aligned}
\label{SomeUEP}
\sum_{l=1}^m\sum_{k\in\Z}\bq_l[n+k]\bq_l[k]=\dde_n=\left\{\begin{array}{cl}
1~&\text{if}~n=0,\vspace{0.4em}\\
0~&\text{if}~n\neq0.
\end{array}\right.\end{aligned}$$ Moreover, once the one dimensional filters generate a wavelet tight frame on $\ell_2(\Z)$, the higher dimensional wavelet tight frame could be obtained via the tensor product of the one dimensional filters.

In the literature, wavelet frames are widely used for the sparse approximation of an image. This is due to the multiscale structure of the wavelet frame systems, short supports of the framelet functions with varied vanishing moments, and the presence of both low pass and high pass filters in the wavelet frame filter banks, which are desirable in sparsely approximating piecewise smooth functions [@B.Dong2017]. For better sparse approximation of an image, the authors in [@J.F.Cai2014] proposed a *data driven tight frame* approach. Specifically, given an image $\bu$, a tight frame system $\bsW$ defined as in \[AnalConv\], which is generated by a set of finitely supported $p\times p$ filters $\{\bq_1,\cdots,\bq_{p^2}\}$ satisfying , is constructed via the following minimization $$\begin{aligned}
\label{DDTFModel}
\min_{\bc,\bsW}~\|\bc-\bsW\bu\|_2^2+\lambda^2\|\bc\|_0~~~~~\text{subject to}~~~~\bsW^T\bsW=\bsI,\end{aligned}$$ with the $\ell_0$ norm $\|\bc\|_0$ encoding the number of nonzero entries in the coefficient vector $\bc$.

Given that each filter is supported on a $p\times p$ grid, we can reformulate \[DDTFModel\] the following way. Reshape all $p\times p$ patches of $\bu$ into $\bU\in\R^{p^2\times m}$ where $m$ denotes the number of total patches. Let $\bsD\in\R^{p^2\times p^2}$ be the matrix generated by concatenating filters $\{\bq_1,\cdots,\bq_{p^2}\}$ into $p^2\times 1$ column vectors $\{\vec{\bq}_1,\cdots,\vec{\bq}_{p^2}\}$. Denote by $\bC\in\R^{p^2\times m}$ the frame coefficients. Hence, we have $$\begin{aligned}
\bu&\Leftrightarrow\bU=(\vec{\bu}_1,\cdots\,\vec{\bu}_m)\in\R^{p^2\times m},\\
\bsW&\Leftrightarrow\bsD=\left(\vec{\bq}_1,\cdots,\vec{\bq}_{p^2}\right)\in\R^{p^2\times p^2},\\
\bc&\Leftrightarrow\bC=(\vec{\bc}_1,\cdots,\vec{\bc}_m)\in\R^{p^2\times m}.\end{aligned}$$ After making these substitutions and relaxing \[SomeUEP\] into $\bsD\bsD^T=p^{-2}\bsI$, \[DDTFModel\] can be reformulated as $$\begin{aligned}
\label{DDTFModelModi}
\min_{\bC,\bsD}~\|\bC-\bsD^T\bU\|_F^2+\lambda^2\|\bC\|_0~~~~~\text{subject to}~~~~\bsD\bsD^T=p^{-2}\bsI,\end{aligned}$$ where $\|\cdot\|_F$ is the Frobenius norm of a matrix. To solve \[DDTFModelModi\], the alternating minimization method with closed form solutions for each stage is presented in [@J.F.Cai2014]. In addition, the proximal alternating minimization (PAM) scheme with global convergence to critical points is proposed in [@C.Bao2015].

Structured low rank matrix approaches in CS-MRI {#SRLACSMRI}
-----------------------------------------------

We begin with introducing some notation. All two dimensional images will be denoted by the bold faced lower case letters and all matrices will be denoted by the bold faced upper case letters. Note that a two dimensional images can also be identified with a vector whenever convenient. We denote by $\OO=\left\{-N/2,\cdots,N/2-1\right\}^2$ with $N\in2\N$ the set of $N\times N$ grid, and the space of complex valued functions on $\OO$ is denoted by $\mI_2\simeq\C^{|\OO|}$. Given two rectangular grids $\KK$ and $\MM$, the contraction $\KK:\MM$ is defined as $$\begin{aligned}
\KK:\MM=\left\{\bk\in\KK:\bk-\bsj\in\KK~\text{for all}~\bsj\in\MM\right\}.\end{aligned}$$ Let $\bsv\in\mI_2$ and let $\KK$ be a rectangular $K_1\times K_2$ grid. The corresponding Hankel matrix $\bmH\bsv$ is an $M_1\times M_2$ matrix ($M_1=|\OO:\KK|$ and $M_2=|\KK|$) generated by concatenating $K_1\times K_2$ patches of $\bsv$ into row vectors. With a slight abuse of notation, we also use $\bmH\w=\left[\begin{array}{cc}
\left(\bmH\w_1\right)^T&\left(\bmH\w_2\right)^T
\end{array}\right]^T$ to denote the $2M_1\times M_2$ two-fold Hankel matrix constructed from $\w=\left(\w_1,\w_2\right)\in\mI_2\times\mI_2$; see \[BlockHankelIllustrate\] for an illustration.

![Illustration of generating a two-fold Hankel matrix $\bmH\w$ from $\w=\left(\w_1,\w_2\right)\in\mI_2\times\mI_2$. Here, $K_1=K_2=2$.[]{data-label="BlockHankelIllustrate"}](BlockHankelIllustrate.pdf){width="90.00000%"}

The conventional CS-MRI approach aims to directly restore the MR image $\bu$ on $\OO$ from the undersampled k-space $\bsf$ satisfying the following two dimensional Fourier transform $$\begin{aligned}
\label{MRForward}
\bsf[\bk]=\wh{u}(L^{-1}\bk)=\int_{\R^2}u(\x)e^{-2\pi i\x\cdot\bk/L}\rd\x,~~~~~~\bk\in\MM,\end{aligned}$$ where $\MM$ is the sampling grid, $L>0$ is the length of field-of-view (FOV), and $u\in L_1(\R^2)$ is the proton spin density distribution in $\R^2$ [@E.M.Haacke1999]. Note that when $\MM=\OO$, i.e. the fully sampled case, $\bu$ is obtained from the inverse DFT $$\begin{aligned}
\label{MRInverse}
\bu[\bp]=\left(\bmsF^{-1}\bsf\right)[\bp]=\f{1}{N^2}\sum_{\bk\in\OO}\bsf[\bk]e^{2\pi i\bp\cdot\bk/N},~~~~~~~\bp\in\OO.\end{aligned}$$ Hence, the conventional CS-MRI approach is based on the DFT $$\begin{aligned}
\label{DFT}
\bsf[\bk]=\left(\bmsF\bu\right)[\bk]=\sum_{\bp\in\OO}\bu[\bp]e^{2\pi i\bp\cdot\bk/N},~~~~~~~\bk\in\MM.
 \end{aligned}$$ Note that this approach is the on-the-grid scheme, as we only consider the information on $\OO$. In other words, we can apply the sparse regularization effectively provided that the singularities of the image $\bu$ are well aligned with $\OO$. However, since the k-space data actually comes from the Fourier transform of a continuous domain function, there exists a basis mismatch between the singularities of $u$ in the continuous domain and the discrete grid $\OO$. Such a basis mismatch would destroy the sparsity structure, leading to a degradation of the restoration performance.

Instead of directly restoring the discrete MR image, the off-the-grid approaches attempt to first restore the fully sampled k-space data $\bsv$ from its undersampled version $\bsf$ modeled as in \[MRForward\]. In other words, we first $$\begin{aligned}
\label{OurProblem}
\text{find}~~\bsv~~~\text{such that}~~\bsv=\msF(u)\big|_{L^{-1}\OO}~~\text{and}~~\bmP_{\MM}\bsv=\bsf,\end{aligned}$$ and then compute $\bu=\bmsF^{-1}\bsv$; see \[CSMRIComparisons\] for the schematic comparison.

![Comparison between the conventional CS-MRI and the off-the-grid CS-MRI approach.[]{data-label="CSMRIComparisons"}](CSMRIComparisons.pdf){width="95.00000%"}

Generally, it is impossible to directly solve \[OurProblem\] without any further information. Nevertheless, by imposing the prior information of $u$ (in the continuous domain) into \[OurProblem\], we can achieve the continuous domain restoration, and reduce the basis mismatch. For this purpose, we follow [@G.Ongie2016] and consider the piecewise constant function model of the proton density $u$: $$\begin{aligned}
\label{uModel}
u(\x)=\sum_{j=1}^J\alpha_j1_{\Om_j}(\x),~~~~~~~\x\in\R^2,\end{aligned}$$ where $\alpha_j\in\C$ and $1_{\Om}$ denotes the characteristic function on a set $\Om$: $1_{\Om}(\x)=1$ if $\x\in\Om$, and $0$ otherwise. Here, we assume that each $\Om_j$ lies in $[-L/2,L/2)^2$. We further assume that \[uModel\] is expressed with the smallest number of characteristic functions such that $\Om_j$’s are pairwise disjoint. Under this setting, the discontinuities of $u$ agrees with $\Gamma=\bigcup_{j=1}^J\p\Om_j$, which is called the *edge set* of $u$ [@G.Ongie2016].

In the sense of distribution, the gradient of $u$ satisfies $$\begin{aligned}
\label{Gradu}
\na u=\left(\p_1 u,\p_2 u\right)=\sum_{j=1}^J\alpha_j\rd\nnu_j,~~~~~~~\text{where}~~~\rd\nnu_j=-\bn_j\rd\ssi\big|_{\p\Om_j},\end{aligned}$$ with $\bn_j$ being the outward normal vector on $\p\Om_j$ and $\ssi$ being the surface measure. Since $|\bn_j|=1$ a.e. on $\p\Om_j$ and $\ssi(\Om_j)<\infty$, $\na u$ in \[Gradu\] defines a finite Radon vector measure on $\R^2$ supported on $\Gamma$. In addition, since \[Gradu\] holds in the sense of tempered distribution, we can compute $\msF(\na u)(\xxi)$ as a Fourier transform of a measure (e.g. [@Folland1999]). Namely, we have $$\begin{aligned}
\label{FouGradu}
\msF(\na u)(\xxi)=\left(\msF(\p_1u)(\xxi),\msF(\p_2u)(\xxi)\right)=\sum_{j=1}^J\alpha_j\int_{\p\Om_j}e^{-2\pi i\xxi\cdot\x}\rd\nnu_j(\x).\end{aligned}$$ Let $\varphi\in\msS$, where $\msS$ is the space of Schwartz functions (e.g. [@Folland1999]). Then the direct computation shows that $\left(\msF(\na u)\ast\wh{\varphi}\right)(\xxi)=\left(\left(\msF(\p_1u)\ast\wh{\varphi}\right)(\xxi),\left(\msF(\p_2u)\ast\wh{\varphi}\right)(\xxi)\right)$ satisfies $$\begin{aligned}
\label{FouGraduConv}
\begin{split}
\left(\msF(\na u)\ast\wh{\varphi}\right)(\xxi)=\sum_{j=1}^J\alpha_j\int_{\p\Om_j}e^{-2\pi i\xxi\cdot\x}\varphi(\x)\rd\nnu_j(\x).
\end{split}\end{aligned}$$ From \[FouGraduConv\], it is easy to see that if $\Gamma\subseteq\left\{\x\in\R^2:\varphi(\x)=0\right\}$, then $\msF(\na u)\ast\wh{\varphi}=0$.

We further assume that there exists a finite (rectangular and symmetric) index set $\KK$ such that $$\begin{aligned}
\label{MajorAssumption}
\Gamma\subseteq\left\{\x\in\R^2:\varphi(\x)=0\right\}~~~~~\text{with}~~\varphi(\x)=\sum_{\bk\in\KK}\a[\bk]e^{2\pi i\bk\cdot\x/L}.\end{aligned}$$ For such $\varphi$, we have $$\begin{aligned}
\label{FouGradConvSemiDiscrete}
\left(\msF(\na u)\ast\wh{\varphi}\right)(\xxi)=\sum_{\bk\in\KK}\msF(\na u)(\xxi-L^{-1}\bk)\a[\bk]=0,~~~~~~~\xxi\in\R^2,\end{aligned}$$ which is referred to as the *linear annihilation relation* [@G.Ongie2015; @G.Ongie2016]. The trigonometric polynomial $\varphi$ in \[MajorAssumption\] is called the *annihilating polynomial*, $\a$ is called the *annihilating filter* with support $\KK$, and $\Gamma$ satisfying \[MajorAssumption\] is called the *trigonometric curves* [@G.Ongie2016]. It is proved in [@G.Ongie2016] that trigonometric curves can always be described as the zero set of a *real valued* trigonometric polynomial, using the tool of algebraic geometries.

In MRI, the Fourier transform of $u$ is sampled on the grid $L^{-1}\OO$ [@E.M.Haacke1999], so \[FouGradConvSemiDiscrete\] becomes the finite system of linear equations $$\begin{aligned}
\label{FouGradConvDiscrete}
\sum_{\bk\in\KK}\msF(\na u)(L^{-1}(\bsl+\bk))\a[-\bk]=0,~~~~~~~\bsl\in\OO:\KK.\end{aligned}$$ In the matrix-vector multiplication form, we have $$\begin{aligned}
\label{RankDeficientHankel}
\bmH\left(\msF(\na u)\big|_{L^{-1}\OO}\right)\a[-\cdot]=\0.\end{aligned}$$ Moreover, if $\varphi$ in \[MajorAssumption\] is the minimal polynomial with the support $\KK$ (i.e. $\KK$ has the smallest linear dimension) and $\KK'$ is the assumed filter support strictly containing $\KK$, then we have $$\begin{aligned}
\label{RankUB}
\rank\left(\bmH\left(\msF(\na u)\big|_{L^{-1}\OO}\right)\right)\leq|\KK'|-|\KK':\KK|.\end{aligned}$$ Roughly speaking, \[RankUB\] means that if $\varphi$ in \[MajorAssumption\] is the minimal polynomial with $\a$ satisfying \[FouGradConvDiscrete,RankDeficientHankel\], so does the translation $\a[\cdot-\bm]$ for all $\bm\in\KK'\setminus\left(\KK':\KK\right)$, or equivalently, $e^{2\pi i\bm\cdot\x}\varphi(\x)$ is also an annihilating polynomial. See [@G.Ongie2015; @G.Ongie2016] for details.

Note that \[RankUB\] means the two-fold Hankel matrix constructed by $\msF(\na u)\big|_{L^{-1}\OO}$ is *rank deficient*. Based on this observation, we can formulate \[OurProblem\] as $$\begin{aligned}
\label{LRSMCompletionRevist}
\min~~\rank\left(\bmH\left(\La\bsv\right)\right)~~~~~~~\text{subject to}~~\bmP_{\MM}\bsv=\bsf.\end{aligned}$$ Here, $\La\in\C^{2|\OO|\times|\OO|}$ is a weight matrix defined as $$\begin{aligned}
\label{Lambda}
\La=\left[\begin{array}{cc}
\La_1&\La_2
\end{array}\right]^T=\left[\begin{array}{cc}
\operatorname{diag}\left(2\pi ik_1/L\right)&\operatorname{diag}\left(2\pi ik_2/L\right)
\end{array}\right]_{\bk=(k_1,k_2)\in\OO}^T,\end{aligned}$$ which is derived from $\msF\left(\na u\right)(\xxi)=2\pi i\xxi\wh{u}(\xxi)$. Hence, the piecewise constant property of $u$ in the continuous domain can be transformed to the low rank property of $\bmH\left(\La\bsv\right)$.

Proposed off-the-grid reconstruction model {#ModelandAlgorithm}
==========================================

Model and algorithm {#ProposedModelandAlgorithm}
-------------------

Denote by $\bsv\in\mI_2$ the k-space data on the grid $\OO$, which is to be restored from a given data $\bsf$ undersampled on $\MM$. Consider a two-fold Hankel matrix $\bmH\left(\La\bsv\right)\in\C^{2M_1\times M_2}$ with $\La$ defined as \[Lambda\]. Note that for $\bq\in\C^{M_2}$, we have $$\begin{aligned}
\bmH\left(\La\bsv\right)\bq=\left(\La\bsv\right)\ast\bq[-\cdot]:=\left[\begin{array}{c}
\left(\La_1\bsv\right)\ast\bq[-\cdot]\vspace{0.25em}\\
\left(\La_2\bsv\right)\ast\bq[-\cdot]
\end{array}\right],\end{aligned}$$ where the two dimensional discrete convolution $\ast$ is performed by reformulating $\bq\in\C^{M_2}$ into a $K_1\times K_2$ filter. Then for an orthogonal matrix $\bY\in\C^{M_2\times M_2}$, we have $$\begin{aligned}
\bmH\left(\La\bsv\right)=\sum_{j=1}^{M_2}\left(\bmH\left(\La\bsv\right)\bY^{(:,j)}\right)\left(\bY^{(:,j)}\right)^*=\sum_{j=1}^{M_2}\left(\left(\La\bsv\right)\ast\bY^{(:,j)}[-\cdot]\right)\left(\bY^{(:,j)}\right)^*.\end{aligned}$$ In addition, it is not hard to see that $$\begin{aligned}
\bmH^*\left(\w\bq^*\right)=\w\ast\overline{\bq}:=\left[\begin{array}{c}
\w_1\ast\overline{\bq}\vspace{0.25em}\\
\w_2\ast\overline{\bq}
\end{array}\right]~~~~~~\text{for}~~~\w=\left[\begin{array}{c}
\w_1\vspace{0.25em}\\
\w_2
\end{array}\right]\in\C^{2M_1}~~\text{and}~~\bq\in\C^{M_2},\end{aligned}$$ and the Moore-Penrose pseudoinverse of $\bmH$ is proportional to $\bmH^*$. Hence, for an orthogonal matrix $\bY\in\C^{M_2\times M_2}$, we have $$\begin{aligned}
\label{HankelTightFrame}
\La\bsv\propto\sum_{j=1}^{M_2}\left(\left(\La\bsv\right)\ast\bY^{(:,j)}[-\cdot]\right)\ast\overline{\bY}^{(:,j)},\end{aligned}$$ where $\propto$ means the pointwise proportionality.

We further introduce $$\begin{aligned}
\bsW&=\left[\bsS_{\bq_1[-\cdot]}^T,\bsS_{\bq_2[-\cdot]}^T,\cdots,\bsS_{\bq_{M_2}[-\cdot]}^T\right]^T,\label{OurAnalysis}\\
\bsW^*&=\left[\bsS_{\overline{\bq}_1},\bsS_{\overline{\bq}_2},\cdots,\bsS_{\overline{\bq}_{M_2}}\right],\label{OurSynthesis}\end{aligned}$$ where $\bq_j=M_2^{-1/2}\bY^{(:,j)}$, and $\bsS_{\bq}$ is defined as in \[DiscreteConv\]. Then up to an appropriate scaling, we can rewrite \[HankelTightFrame\] as $$\begin{aligned}
\label{OrthogonalFilterBanks}
\La\bsv=\sum_{j=1}^{M_2}\bsS_{\overline{\bq}_j}\left(\bsS_{\bq_j[-\cdot]}\left(\La\bsv\right)\right)=\bsW^*\bsW\left(\La\bsv\right),\end{aligned}$$ which implies that the orthogonal matrix $\bY\in\C^{M_2\times M_2}$ indeed generates a set of filter banks of a tight frame system.

To reconstruct $\bsv=\msF(u)\big|_{L^{-1}\OO}$ corresponding to a piecewise constant function $u$ in \[uModel\], we assume that $\rank\left(\bmH\left(\La\bsv\right)\right)=r\ll2M_1\wedge M_2$, following \[SRLACSMRI\]. Consider its SVD $$\begin{aligned}
\label{HankelSVD}
\bmH\left(\La\bsv\right)=\left[\begin{array}{c}
\bmH\left(\La_1\bsv\right)\\
\bmH\left(\La_2\bsv\right)
\end{array}\right]=\bX\Sig\bY^*=\sum_{j=1}^r\Sig^{(j,j)}\bX^{(:,j)}\left(\bY^{(:,j)}\right)^*,\end{aligned}$$ with $\Sig^{(1,1)}\geq\Sig^{(2,2)}\geq\cdots\geq\Sig^{(r,r)}>0$ and $\Sig^{(j,j)}=0$ for $j>r$. Since we have $$\begin{aligned}
\label{HankelSVDConv}
\Sig^{(j,j)}\bX^{(:,j)}=\bmH\left(\La\bsv\right)\bY^{(:,j)}=\left(\La\bsv\right)\ast\bY^{(:,j)}[-\cdot],\end{aligned}$$ it follows that $$\begin{aligned}
\label{HankelSVDSparse}
\begin{split}
\left(\La\bsv\right)\ast\bY^{(:,j)}[-\cdot]&\neq\0,~~~~~~~j=1,\cdots,r,\vspace{0.25em}\\
\left(\La\bsv\right)\ast\bY^{(:,j)}[-\cdot]&=\0,~~~~~~~j=r+1,\cdots,M_2,
\end{split}\end{aligned}$$ which leads to $$\begin{aligned}
\label{LowRankHankelSVD}
\La\bsv\propto\sum_{j=1}^r\left(\left(\La\bsv\right)\ast\bY^{(:,j)}[-\cdot]\right)\ast\overline{\bY}^{(:,j)}.\end{aligned}$$ Hence, if we define $\bsW$ and $\bsW^*$ as in \[OurAnalysis,OurSynthesis\] using $\bY\in\C^{M_2\times M_2}$ in \[HankelSVD\], we have $$\begin{aligned}
\label{TightFrameSparse}
\begin{split}
\bsS_{\bq_j[-\cdot]}\left(\La\bsv\right)&\neq\0,~~~~~~~j=1,\cdots,r,\vspace{0.25em}\\
\bsS_{\bq_j[-\cdot]}\left(\La\bsv\right)&=\0,~~~~~~~j=r+1,\cdots,M_2,
\end{split}\end{aligned}$$ which finally leads to $$\begin{aligned}
\label{LowRankHankelTightFrame}
\La\bsv=\sum_{j=1}^r\bsS_{\overline{\bq}_j}\left(\bsS_{\bq_j[-\cdot]}\left(\La\bsv\right)\right)=\bsW^*\bsW\left(\La\bsv\right).\end{aligned}$$ This means that, if $\rank\left(\bmH\left(\La\bsv\right)\right)=r$, then its SVD constructs a tight frame such that $\La\bsv$ has only $r$ *nonzero* canonical coefficients.

Based on these observations, we can consider constructing an *adaptive* tight frame under which $\La\bsv$ has a small number of canonical nonzero coefficients, instead of directly minimizing $\rank\left(\bmH\left(\La\bsv\right)\right)$. However, it is in general difficult to directly impose such an idea on the reconstruction model for the following reasons: 1) it is an NP-hard problem to minimize the number of nonzero coefficients that represents/approximates $\La\bsv$ even under a given tight frame; 2) the penalization on canonical coefficients in general requires an iterative solver, which will be too expensive in constructing an adaptive tight frame [@J.F.Cai2014]; 3) the available k-space data is always corrupted by the noise in practice, which makes \[LowRankHankelTightFrame\] only approximately true in practice. Nevertheless, \[LowRankHankelTightFrame\] implies that there exists an adaptive tight frame under which $\La\bsv$ can be *sparsely approximated*, so we can adopt the idea of data driven tight frames [@J.F.Cai2014]. In particular, we propose the following data driven tight frame (DDTF) based continuous domain CS-MRI reconstruction model $$\begin{aligned}
\label{ProposedCSMRIModel}
\begin{split}
&~~~\min_{\bsv,\bc,\bsW}~\f{1}{2}\left\|\bmP_{\MM}\bsv-\bsf\right\|_2^2+\f{\mu}{2}\left\|\bsW\left(\La\bsv\right)-\bc\right\|_2^2+\gamma\left\|\bc\right\|_0\\
&\text{subject to}~~\bsv\in\mC~~\text{and}~~\bsW^*\bsW=\bsI
\end{split}\end{aligned}$$ to reconstruct the fully sampled k-space data $\bsv$ and then obtain the MR image by $\bu=\bmsF^{-1}\bsv$. Here, $\bsW$ is a tight frame transform defined as in \[OurAnalysis\] whose filter banks have to be learned, and $\mC$ is a constraint set imposing the boundedness of $\bsv$. In this paper, we choose $$\begin{aligned}
\label{ConstraintSet}
\mC=\left\{\bsv\in\mI_2:|\bsv[\bk]|\leq R~~~\text{for all}~~\bk\in\OO\right\}\simeq\DD_R^{N^2},\end{aligned}$$ where $\DD_{R}=\left\{v\in\C:|v|\leq R\right\}$ with a sufficiently large $R>0$.

\[RK1\] The choice of the constraint set $\mC$ comes from the following: if $u\in L_1(\R^2)$ is modeled as in \[uModel\], then from \[FouGradu\], we have $$\begin{aligned}
\left|\msF(\na u)(\xxi)\right|\leq\sum_{j=1}^J|\alpha_j|\ssi(\Om_j)<\infty~~~~~\text{for all}~~\xxi\in\R^2.\end{aligned}$$ In other words, $\wh{u}(\xxi)$ decays at infinity in the sense that $$\begin{aligned}
\sup_{\xxi\in\R^2}\left|\xxi\right|\left|\wh{u}(\xxi)\right|<\infty,\end{aligned}$$ and this implies that $\wh{u}(\xxi)$ (and whence $\bsv[\bk]$ since we want $\bsv=\msF(u)\big|_{L^{-1}\OO}$) has to be bounded. In addition, as we shall see in \[ConvergenceAnalysis\], the boundedness constraint is required for the convergence analysis of the proximal alternating scheme to solve \[ProposedCSMRIModel\]. That being said, numerically, the constraint set has a very minor effect on the restoration results provided that $R>0$ is sufficiently large. In fact, the proximal alternating minimization seems to converge even without using $\mC$.

In \[ProposedCSMRIModel\], the sparsity promoting $\ell_0$ norm $\left\|\bc\right\|_0$ is used as a relaxation of \[TightFrameSparse\] (or equivalently, \[HankelSVDSparse\]), and the tight frame constraint $\bsW^*\bsW=\bsI$ is related to \[LowRankHankelTightFrame\]. Finally, the term $\left\|\bsW\left(\La\bsv\right)-\bc\right\|_2^2$ is to provide the flexibility in sparse approximation. Namely, given that $\bsW$ is a tight frame, we have $$\begin{aligned}
\left\|\bsW\left(\La\bsv\right)-\bc\right\|_2^2=\left\|\bsW^*\bc-\La\bsv\right\|_2^2+\left\|\left(\bsI-\bsW\bsW^*\right)\bc\right\|_2^2,\end{aligned}$$ which yields a balanced approach (e.g. [@J.F.Cai2008; @R.H.Chan2003]). In this work, we choose this balanced approach for the following two reasons. Firstly, the k-space data is in general corrupted by the (thermal) noise [@E.M.Haacke1999] and $\La$ amplifies the noise in the high frequency, which always hampers $\bsW\left(\La\bsv\right)$ from being sparse (or $\bmH\left(\La\bsv\right)$ being low rank). Secondly, the direct penalization on the canonical coefficient $\bsW\left(\La\bsv\right)$ requires an iterative solver, which will be too expensive in constructing an adaptive tight frame [@J.F.Cai2014]. Therefore, by promoting the distance between $\bsW\left(\La\bsv\right)$ and $\bc$, we expect to achieve a better sparse approximation of $\La\bsv$ (and thus a better low rank approximation of $\bmH\left(\La\bsv\right)$) with a low computational complexity.

To solve \[ProposedCSMRIModel\], we use the proximal alternating minimization (PAM) algorithm introduced in [@H.Attouch2010]. We initialize $\bsv_0=\bmP_{\mC}\left(\bmP_{\MM}\bsf\right)$ where $$\begin{aligned}
\label{ProjectOntoC}
\bmP_{\mC}\left(\bsv\right)[\bk]=\min\left\{\left|\bsv[\bk]\right|,R\right\}\exp\left\{i\arg\bsv[\bk]\right\}.\end{aligned}$$ The initializations of $\bsW_0$ and $\bc_0$ are obtained from the SVD of $\bmH\left(\La\bsv_0\right)=\bX_0\Sig_0\bY_0^*$: $$\begin{aligned}
\label{Initialize}
\bsW_0=M_2^{-1/2}\bY_0~~~~~\text{and}~~~~~\bc_{0,j}=\left\{\begin{array}{cl}
M_2^{-1/2}\left(\La\bsv_0\right)\ast\bY_0^{(:,j)}[-\cdot]~&j=1,\cdots,r,\vspace{0.4em}\\
\0~&j=r+1,\cdots,M_2,
\end{array}\right.\end{aligned}$$ where $r$ is an estimated rank. After the initializations, we optimize $(\bsv,\bc,\bsW)$ alternatively, as summarized in \[Alg1\].

It is easy to see that each subproblem in \[Alg1\] has a closed form solution. The solution to \[vsubprob\] is given by $$\begin{aligned}
\label{vsubexplicit}
\begin{split}
\bsv_{n+1/2}&=\left(\bmP_{\MM}+\mu\La^*\La+\beta_1\bsI\right)^{-1}\left(\bmP_{\MM}\bsf+\mu\La^*\bsW_n^*\bc_n+\beta_1\bsv_n\right),\\
\bsv_{n+1}&=\bmP_{\mC}\left(\bsv_{n+1/2}\right),
\end{split}\end{aligned}$$ where $\bmP_{\mC}$ is defined as \[ProjectOntoC\]. It is worth noting that since $\bmP_{\MM}+\mu\La^*\La+\beta_1\bsI$ is a diagonal matrix, no matrix inversion is needed.

To solve \[Wsubprob,csubprob\], we introduce the following substitutions $$\begin{aligned}
\label{PatchReformulation}
\left(\La\bsv,\bc,\bsW\right)\Leftrightarrow\left(\bsH,\bC,\bsD\right),\end{aligned}$$ where $\bsH=\left[\begin{array}{cc}
\bsH_1&\bsH_2
\end{array}\right]$ with each $\bsH_j$ for $j=1$ and $2$ generated by using $K_1\times K_2$ patches of $\La_j\bsv$. Then \[Wsubprob,csubprob\] can be reformulated as $$\begin{aligned}
\bC_{n+1}&={\operatornamewithlimits{argmin}}~\gamma\left\|\bC\right\|_0+\f{\mu}{2}\left\|\bC-\bsD_{n}^T\bsH_{n+1}\right\|_F^2+\f{\beta_2}{2}\left\|\bC-\bC_n\right\|_F^2,\label{Cupdate}\\
\bsD_{n+1}&={\operatornamewithlimits{argmin}}_{\bsD\bsD^*=M_2^{-1}\bsI}~\f{\mu}{2}\left\|\bsD^T\bsH_{n+1}-\bC_{n+1}\right\|_F^2+\f{\beta_3}{2}\left\|\bsD-\bsD_n\right\|_F^2,\label{Dupdate}\end{aligned}$$ respectively. The closed form solution to \[Cupdate\] is expressed as $$\begin{aligned}
\label{csubexplicit}
\bC_{n+1}=\bmT_{\sqrt{2\gamma/(\mu+\beta_2)}}\left(\f{\mu\bsD_{n}^T\bsH_{n+1}+\beta_2\bC_n}{\mu+\beta_2}\right),\end{aligned}$$ where $\bmT_{\gamma}$ is a hard thresholding; $\bmT_{\gamma}\left(\bC\right)[m,n]=\bC[m,n]$ if $|\bC[m,n]|>\gamma$, and $0$ otherwise. To solve \[Dupdate\], we use the following closed form formula $$\begin{aligned}
\label{Wsubexplicit}
\bsD_{n+1}=M_2^{-1/2}\bX_n\bY_n^*~~~\text{where}~~\bX_n\Sig_n\bY_n^*=\overline{\bsH}_{n+1}\bC_{n+1}^T+\f{\beta_3}{\mu}\bsD_n~~\text{is the SVD}.\end{aligned}$$ Hence, we only compute the SVD of an $M_2\times M_2$ matrix $\overline{\bsH}_{n+1}\bC_{n+1}^T+\left(\beta_3/\mu\right)\bsD_n$, leading to the computational efficiency over directly minimizing the rank of a $2M_1\times M_2$ matrix $\bmH\left(\La\bsv\right)$.

We further mention that it is not necessary to explicitly perform \[PatchReformulation\]. Noting that $\bsH$ is the transpose of $\bmH\left(\La\bsv\right)$, we can compute $\bsD_{n}^T\bsH_{n+1}$ in \[csubexplicit\] and $\overline{\bsH}_{n+1}\bC_{n+1}^T$ in \[Wsubexplicit\] by using $3M_2$ ($M_2$ for $\bsD_{n}^T\bsH_{n+1}$ and $2M_2$ for $\overline{\bsH}_{n+1}\bC_{n+1}^T$) fast Hankel matrix-vector multiplications/two dimensional convolutions directly from $\La\bsv_{n+1}$, requiring $O(M_2N^2\log N)$ floating point operations. In addition, since $\bmT_{\gamma}$ is an elementwise operator, we can update each column of $\overline{\bsH}_{n+1}\bC_{n+1}^T+\left(\beta_3/\mu\right)\bsD_n$ in \[Wsubexplicit\] directly after updating each row of $\bC_{n+1}$ by \[csubexplicit\].

Convergence analysis {#ConvergenceAnalysis}
--------------------

In this subsection, we focus on the convergence analysis of the sequence $\left(\bsv_n,\bc_n,\bsW_n\right)$ generated by \[Alg1\] based on the frameworks in [@H.Attouch2010; @C.Bao2016; @J.Bolte2014; @Kurdyk1998; @Lojasiewicz1995; @Y.Xu2013] for the nonconvex and nondifferentiable optimization. To do this, we begin with introducing some basic notation and definitions.

Let $f:\R^n\to\R\cup\{\infty\}$ be a proper and lower semicontinuous (lsc) function.

**

1.  The domain of $f$, denoted $\dom(f)$, is defined as $$\begin{aligned}
    \dom(f)=\{\bu\in\R^n:f(\bu)<\infty\}.\end{aligned}$$

2.  For each $\bu\in\dom(f)$, the Fréchet subdifferential of $f$ at $\bu$ is defined as $$\begin{aligned}
    \p_Ff(\bu)=\left\{\bs\in\R^n:\liminf_{\bsv\to\bu}\f{f(\bsv)-f(\bu)-\la\bs,\bsv-\bu\ra}{\|\bsv-\bu\|}\geq0\right\}.\end{aligned}$$ If $\bu\notin\dom(f)$, then we set $\p_Ff(\bu)=\emptyset$.

3.  The (limiting-) subdifferential of $f$ at $\bu$ is defined as $$\begin{aligned}
    \p f(\bu)=\left\{\bs\in\R^n:\exists \bu_n~\text{s.t.}~f(\bu_n)\to f(\bu)~\&~\bs_n\in\p_Ff(\bu_n)\to\bs\right\}.\end{aligned}$$ The domain of $\p f$ is defined as $\dom(\p f)=\left\{\bu\in\R^n:\p f(\bu)\neq\emptyset\right\}$.

4.  $\bu\in\dom(f)$ is a critical point of $f$ if $\0\in\p f(\bu)$.

\[RK3\] Since the minimization of a real-valued complex variable objective function in \[ProposedCSMRIModel\] is based on the Wirtinger calculus (i.e. by identifying $\C\simeq\R^2$), the theoretical results on the proximal schemes for the optimization on the real vector space are still applicable.

The global convergence of \[Alg1\] is based on the Kurdyka-Lojasiewicz (KL) property in [@Kurdyk1998; @Lojasiewicz1995]. It is in general difficult to directly verify whether a given function $f$ satisfies the KL property by the definition. Instead, we can verify this through several special types of functions which are proven to satisfy the KL property, such as analytic functions and semi-algebraic functions (e.g. [@H.Attouch2010; @C.Bao2016; @J.Bolte2014; @Y.Xu2013]).

For the notational simplicity, we introduce $\z=\left(\bsv,\bc,\bsW\right)$ and $$\begin{aligned}
P(\z)=\f{1}{2}\left\|\bmP_{\MM}\bsv-\bsf\right\|_2^2+\f{\mu}{2}\left\|\bsW\left(\La\bsv\right)-\bc\right\|_2^2.\end{aligned}$$ Recall that $\mC\simeq\DD_{R}^{N^2}$ as in \[ConstraintSet\], and we introduce $$\begin{aligned}
\mD=\left\{\bsW\in\C^{M_2N^2\times N^2}:\bsW^*\bsW=\bsI_{N^2\times N^2}\right\}.\end{aligned}$$ Using this notation, we reformulate \[ProposedCSMRIModel\] as $$\begin{aligned}
\label{ObjectModi}
\min_{\z}~H(\z):=P(\z)+\gamma\left\|\bc\right\|_0+\imath_{\mC}(\bsv)+\imath_{\mD}(\bsW),\end{aligned}$$ where $\imath_{\mA}$ is the indicator function of a set $\mA$: $\imath_{\mA}(\z)=0$ if $\z\in\mA$, and $\infty$ otherwise.

It is not hard to verify that the indicator functions are analytic, and $\left\|\cdot\right\|_0$ is a semi-algebraic function. Since the remaining term $P(\z)$ is merely a polynomial, we can conclude that the function $H(\z)$ in \[ObjectModi\] satisfies the KL-property. Under this formulation, we can present the following result on the global convergence of \[Alg1\].

\[Th1\] Let $H(\z)$ be defined as \[ObjectModi\]. Then the sequence $\{\z_n:n\in\N\}$ generated by \[Alg1\] is globally convergent, and its limit is a critical point of $H$. Moreover, we have $$\begin{aligned}
\label{FiniteLength}
\sum_{n=0}^{\infty}\left\|\z_{n+1}-\z_n\right\|_2<\infty,\end{aligned}$$ i.e. the sequence $\{\z_n:n\in\N\}$ satisfies the finite length property.

Since the proof of \[Th1\] follows the framework given in [@C.Bao2016 Theorem 3.7.], we only need to mention the key component. Since $H(\z)\geq0$ satisfies the KL property, it can be easily verified that it satisfies the first condition in [@C.Bao2016 Theorem 3.7]. Interested readers may refer to [@C.Bao2016] for the details. Hence, we can complete the proof provided that we can verify that the sequence $\left\{\z_n:n\in\N\right\}$ generated by \[Alg1\] is bounded. In fact, since the constraints $\mC$ and $\mD$ are compact sets, it suffices to show that $\bc_n$ is bounded.

\[Lemma1\] Let $H(\z)$ be defined as in . For the sequence $\left\{\z_n:n\in\N\right\}$ generated by \[Alg1\], there exists a constant $C>0$ independent on $n$ such that $$\begin{aligned}
\left\|\bc_n\right\|_2\leq C~~~~~~~\text{for all}~~~n\in\N.\end{aligned}$$ Hence, $\bc_n$ is bounded.

Since $\bsv_n\in\mC$ for all $n\geq0$, we have $\left\|\bsv_n\right\|_2\leq RN$ for all $n\geq0$. Hence, the proof is completed by mathematical induction. For $n=0$, \[Initialize\] implies $$\begin{aligned}
\left\|\bc_0\right\|_2\leq\left\|\bsW_0\left(\La\bsv_0\right)\right\|_2=\left\|\La\bsv_0\right\|_2\leq\left\|\La\right\|\left\|\bsv_0\right\|_2\leq\left\|\La\right\|RN,\end{aligned}$$ where $\left\|\La\right\|$ denotes the spectral norm of $\La$, and the equality comes from the fact that $\bsW_0\in\mD$.

For the mathematical induction, we assume that $\left\|\bc_n\right\|_2\leq C:=\left\|\La\right\|RN$. Due to the hard thresholding in \[csubexplicit\], we have $$\begin{aligned}
\left\|\bc_{n+1}\right\|_2&\leq\left\|\f{\mu\bsW_{n}\left(\La\bsv_{n+1}\right)+\beta_2\bc_n}{\mu+\beta_2}\right\|_2\\
&\leq\f{\mu}{\mu+\beta_2}\left\|\bsW_{n}\left(\La\bsv_{n+1}\right)\right\|_2+\f{\beta_2}{\mu+\beta_2}\left\|\bc_n\right\|_2\\
&=\f{\mu}{\mu+\beta_2}\left\|\La\bsv_{n+1}\right\|_2+\f{\beta_2}{\mu+\beta_2}\left\|\bc_n\right\|_2\\
&\leq\f{\mu}{\mu+\beta_2}\left\|\La\right\|\left\|\bsv_{n+1}\right\|_2+\f{\beta_2}{\mu+\beta_2}\left\|\bc_n\right\|_2\leq\left\|\La\right\|RN=C,\end{aligned}$$ where the equality again comes from the fact that $\bsW_n\in\mD$ for all $n\in\N$. This completes the proof.

Experimental results {#Experiments}
====================

In this section, we present the experimental results on the phantom image and the real MR image used in [@G.Ongie2016], to compare the proposed DDTF based CS-MRI model \[ProposedCSMRIModel\] with several existing methods. Note that the major focus of this paper is to reconstruct a piecewise constant MR image from a given undersampled k-space data. For this purpose, we choose to compare with the total variation (TV) model [@M.Lustig2007] $$\begin{aligned}
\label{TVModel}
\min_{\bu}~\f{1}{2}\left\|\bmP_{\MM}\bmsF\bu-\bsf\right\|_2^2+\gamma\left\|\Na\bu\right\|_1,\end{aligned}$$ and the Haar framelet (Haar) model (e.g. [@J.F.Cai2012]) $$\begin{aligned}
\label{HaarModel}
\min_{\bu}~\f{1}{2}\left\|\bmP_{\MM}\bmsF\bu-\bsf\right\|_2^2+\left\|\gga\cdot\bsW\bu\right\|_1,\end{aligned}$$ where $\gga=\left\{\gamma 2^{-l}:l=0,\cdots,L-1\right\}$. Both \[TVModel,HaarModel\] are solved by the split Bregman method [@J.F.Cai2009/10; @T.Goldstein2009]. We also compare with the following Schatten $p$-norm minimization model $$\begin{aligned}
\label{Schatten}
\min_{\bsv\in\mC}~\left\|\bmP_{\MM}\bsv-\bsf\right\|_2^2+\gamma\left\|\bmH\left(\La\bsv\right)\right\|_p^p,\end{aligned}$$ with $p=0$, $0.5$, and $1$, solved by the generic iterative reweighted annihilating filters (GIRAF) method [@G.Ongie2017] based on the split Bregman algorithm, which are referred to as “GIRAF$0$”, “GIRAF$0.5$”, and “GIRAF$1$”, respectively. All experiments are implemented on MATLAB $\mathrm{R}2014\mathrm{a}$ running on a laptop with $64\mathrm{GB}$ RAM and Intel(R) Core(TM) CPU $\mathrm{i}7$-$8750\mathrm{H}$ at $2.20\mathrm{GHz}$ with $6$ cores.

      Dataset                 Model                $K$       $r$     $\mu$    $\gamma$    $\beta_1$   $\beta_2$   $\beta_3$  
  --------------- ----------------------------- --------- --------- -------- ----------- ----------- ----------- ----------- --
                         TV \[TVModel\]          $\cdot$   $\cdot$    $10$     $0.05$      $\cdot$     $\cdot$     $\cdot$   
                       Haar \[HaarModel\]        $\cdot$   $\cdot$    $10$     $0.025$     $\cdot$     $\cdot$     $\cdot$   
                       GIRAF \[Schatten\]         $25$     $\cdot$    $4$     $10^{-6}$    $\cdot$     $\cdot$     $\cdot$   
                   DDTF \[ProposedCSMRIModel\]    $25$      $500$    $0.1$      $10$      $10^{-4}$   $10^{-4}$   $10^{-4}$  
   \*[Real MR]{}         TV \[TVModel\]          $\cdot$   $\cdot$    $10$     $0.01$      $\cdot$     $\cdot$     $\cdot$   
                       Haar \[HaarModel\]        $\cdot$   $\cdot$    $10$     $0.025$     $\cdot$     $\cdot$     $\cdot$   
                       GIRAF \[Schatten\]         $33$     $\cdot$    $2$     $10^{-6}$    $\cdot$     $\cdot$     $\cdot$   
                   DDTF \[ProposedCSMRIModel\]    $33$      $871$    $0.05$      $5$      $10^{-4}$   $10^{-4}$   $10^{-4}$  

  : Parameter selection for each dataset. For TV \[TVModel\], Haar \[HaarModel\], and GIRAF \[Schatten\], $\mu$ refers to the internal parameter for the split Bregman algorithm.[]{data-label="TableParameter"}

In all cases, we choose $R=|\bsf[\0]|$ if $\0\in\MM$ and $R=10^8$ otherwise for the constraint set $\mC$ in \[ConstraintSet\]. We use the forward difference for the discrete gradient $\Na$ in \[TVModel\], and $\bsW$ in \[HaarModel\] is chosen to be the undecimated tensor product Haar framelet transform with $1$ level of decomposition [@B.Dong2013]. Both \[ProposedCSMRIModel,Schatten\] use the $K\times K$ square patch to generate the two-fold Hankel matrix for simplicity. As a rule of a thumb, we choose $K$ to be around $9\sim13\%$ of $N$, and we choose $r\approx0.8K^2$ for the initialization \[Initialize\] of \[ProposedCSMRIModel\], both of which depend on the geometry of the target image. The detailed choice of the remaining regularization parameters are summarized in \[TableParameter\]. Empirically, we observe that $\mu\approx0.01\gamma$ is an appropriate choice for \[ProposedCSMRIModel\], and we further observe that when $\mu$ is large, the restored k-space data $\bsv$ has a faster decay than smaller $\mu$. Hence, the parameters are manually tuned so that we can achieve the optimal restoration of both the low frequencies and high frequencies. For the on-the-grid approaches \[TVModel,HaarModel\], the stopping criterion is $$\begin{aligned}
\label{StoppingOntheGrid}
\f{\left\|\bu_{n+1}-\bu_n\right\|_2}{\left\|\bu_n\right\|_2}\leq{\varepsilon}:=2\times10^{-4},\end{aligned}$$ and for the off-the-grid approaches \[ProposedCSMRIModel,Schatten\], we use $$\begin{aligned}
\label{StoppingOfftheGrid}
\f{\left\|\bsv_{n+1}-\bsv_n\right\|_2}{\left\|\bsv_n\right\|_2}\leq{\varepsilon}.\end{aligned}$$ We also set the maximum allowable number of iterations to be $600$. To measure the quality of restored images, we compute the signal-to-noise ratio (SNR) and the high frequency error norm (HFEN) [@S.Ravishankar2011] of the restored images. Note that for the off-the-grid approaches \[ProposedCSMRIModel,Schatten\], the restored image is computed via the inverse DFT of the restored k-space data; see \[CSMRIComparisons\].

Phantom experiments {#PhantomExperiments}
-------------------

For the piecewise constant phantom experiments, we first compute the fully sampled k-space data from the analytical frequency domain expressions of the MR phantoms, as performed in [@M.Guerquin-Kern2012]. Then using the variable density random sampling method in [@M.Lustig2007], we generate $20\%$ undersampled k-space data. The complex white Gaussian noise is also added so that the resulting SNR of the samples is approximately $25\mathrm{dB}$ (See \[PhantomDataSet\]).

  ------ ---------- ---------- ---------- ---------- ---------- ---------- -------------------
                                                                           
                                            $p=0$     $p=0.5$     $p=1$    
   SNR     $9.23$    $20.29$    $20.99$    $21.12$    $20.88$    $17.13$    $\textbf{26.69}$
   HFEN   $0.5490$   $0.0992$   $0.0904$   $0.0948$   $0.0999$   $0.1828$   $\textbf{0.0572}$
  ------ ---------- ---------- ---------- ---------- ---------- ---------- -------------------

  : Comparison of signal-to-noise ratio and high frequency error norm for the phantom experiments.[]{data-label="PhantomTable"}

\[PhantomTable\] summarizes the SNR and the HFEN of the aforementioned restoration models, and \[PhantomResults\] displays the visual comparisons with the zoom-in views in \[PhantomResultsZoom\] and the error maps in \[PhantomResultsErrorMap\], respectively. We can see that the proposed data driven tight frame model \[ProposedCSMRIModel\] consistently outperforms both the on-the-grid approaches (\[TVModel,HaarModel\]) and the existing off-the-grid approaches \[Schatten\] with a smaller error map in \[PhantomDDTFError\]. Noting that \[ProposedCSMRIModel\] is an off-the-grid approach, the experimental results also suggest that the off-the-grid approaches have better performance in the CS-MRI due to its ability to reducing the basis mismatch between the true support (or the true singularity) in continuum and the discrete grid. In fact, due to this basis mismatch, we can see from \[PhantomTV,PhantomHaar,PhantomTVZoom,PhantomHaarZoom\] that the on-the-grid approaches lead to the distortions of three small ellipses, and the errors concentrate on the edges (\[PhantomTVError,PhantomHaarError\]) compared to the off-the-grid approaches.

It is also worth noting that among the off-the-grid approaches, the proposed DDTF model introduces less artifacts near the edges. In the literature, the noise in the k-space data affects $\rank\left(\bmH\left(\La\bsv\right)\right)$ even in the fully sampled case as the weight matrix $\La$ amplifies the noise in the high frequencies. Hence under such an amplified noise, it is likely that the direct rank minimization leads to the artifacts near the edges corresponding to the high frequencies in the frequency domain, as shown in \[PhantomGIRAF0,PhantomGIRAFHalf,PhantomGIRAF1,PhantomGIRAF0Zoom,PhantomGIRAFHalfZoom,PhantomGIRAF1Zoom,PhantomGIRAF0Error,PhantomGIRAFHalfError,PhantomGIRAF1Error\]. In contrast, the sparse approximation of $\La\bsv$ can achieve the denoising effect in spite of the amplified noise, leading to the better restoration results with less artifacts near the edges.

For further comparisons, we also present the restored k-space data (in the log scale) in \[PhantomResultsk\]. Note that since the sampling is dense in the low frequencies while the high frequencies are loosely sampled, the restoration qualities depend heavily on the restoration accuracy of high frequency k-space data. Indeed, we can see from \[PhantomTVk,PhantomHaark\] that the restored k-space data by \[TVModel,HaarModel\] decays faster than the original one, which also leads to the inferior restoration results. Even though the GIRAF models are in general able to restore the high frequency part better than the on-the-grid approaches, they still fail to restore the dominant structures on the high frequencies, as shown in \[PhantomGIRAF0k,PhantomGIRAFHalfk,PhantomGIRAF1k\]. In contrast, the proposed model \[ProposedCSMRIModel\] is able to restore the high frequency k-space data in spite of the loose sampling, which also results in the improvements over the existing approaches. In summary, our proposed DDTF CS-MRI model shows the overall better restoration quality in both the indices (SNR and HFEN) and the visual quality.

\

\

\

\

Real MR image experiments {#RealMRExperiments}
-------------------------

The real MR image experiments use the k-space data which is obtained from a fully sampled $4$-coil acquisition, and then compressed into a single virtual coil using the SVD technique in [@T.Zhang2013]. Since the data from the single virtual coil is complex-valued in the image domain with smoothly varying phase, we further correct the phase using the method described in [@G.Ongie2016]. More concretely, we first perform the inverse DFT of the zero padded k-space data, canceling out the phase in the image domain, and passing back to the frequency domain. Then as in the phantom experiments, we generate $20\%$ undersampled k-space data using the variable density sampling, and further add the complex white Gaussian noise so that the resulting SNR of the samples is approximately $25\mathrm{dB}$; see \[RealDataSet\].

  ------ ---------- ---------- ---------- ---------- ---------- ---------- -------------------
                                                                           
                                            $p=0$     $p=0.5$     $p=1$    
   SNR    $11.09$    $14.27$    $14.31$    $15.36$    $15.21$    $13.25$    $\textbf{16.15}$
   HFEN   $0.6050$   $0.3477$   $0.3393$   $0.2924$   $0.3049$   $0.4409$   $\textbf{0.2521}$
  ------ ---------- ---------- ---------- ---------- ---------- ---------- -------------------

  : Comparison of signal-to-noise ratio and high frequency error norm for the real MR experiments.[]{data-label="RealTable"}

\[RealTable\] summarizes the SNR and the HFEN of the restoration results. For visual comparisons, the restored images, the error maps, and the restored k-space data are presented in \[RealResults,RealResultsk,RealResultsErrorMap\], respectively. Overall, we can see that the pros and the cons of the restoration methods are similar to the phantom experiments, and our proposed model \[ProposedCSMRIModel\] again consistently outperforms the other restoration methods. It is also worth noting that the proposed DDTF model restores more structured k-space data than the other models. In fact, the real k-space data is in general contaminated by the thermal noise [@E.M.Haacke1999], which makes the fully sampled k-space data less structured than the piecewise constant phantom. Together with the noise amplification by $\La$ on the high frequency part, even the fully sampled k-space data may not correspond to a low rank two-fold Hankel matrix. Nevertheless, we can observe that despite such noise amplifications, the sparse approximation by the data driven tight frame is likely to achieve better denoising than other models, leading to the better restoration results. Finally, we list some zoom-in views in \[RealResultsZoom\] to illustrate that our model \[ProposedCSMRIModel\] can restore structures better than the other existing methods.

\

\

\

\

Conclusion and future directions {#Conclusion}
================================

In this paper, we propose a new off-the-grid CS-MRI reconstruction model for the piecewise constant image restoration in the two dimensional FRI framework [@G.Ongie2016]. Our proposed model is inspired by the observation that the SVD of a Hankel matrix to some extent corresponds to an adaptive tight frame system which can represent a given image with a small number of nonzero canonical coefficients. This motivates us to adopt the sparse approximation by the data driven tight frames as an alternative to the structured low rank matrix completion. Finally, the numerical experiments show that our approach outperforms both the conventional on-the-grid approaches and the existing the low rank Hankel matrix models and their relaxations. To solve the nonconvex and nondifferentiable model, a proximal alternating minimization algorithm is presented, which has guaranteed global convergence to a critical point.

For the future work, we plan to provide a rigorous theoretical framework on our observations. Specifically, we need to rigorously analyze whether the data driven tight frame can indeed reflect the true frequency information of spectrally sparse signals (Such signals correspond to low rank Hankel matrices). In addition, it is also likely to extend the idea in this paper to the piecewise smooth image restoration framework, such as the total generalized variation [@K.Bredies2010] and the combined first and second order TV model [@M.Bergounioux2010; @K.Papafitsoros2014], by considering the higher order derivatives. It would be also interesting to find other pseudodifferential operators/Fourier integral operators which can provide more insightful information on the rank of the structured matrix constructed from the Fourier transform of a piecewise constant function. For example, we can attempt to find transformations under which the rank of the structured matrix is related to $J$ in \[uModel\].

Acknowledgments {#acknowledgments .unnumbered}
===============

The authors would like to thank Dr. Greg Ongie in the Department of Statistics at the University of Chicago, an author of [@G.Ongie2018; @G.Ongie2015; @G.Ongie2016; @G.Ongie2017], for making the data sets as well as the MATLAB toolbox available so that the experiments can be implemented.

[10]{}

, [*Proximal alternating minimization and projection methods for nonconvex problems: an approach based on the [K]{}urdyka-[L]{}ojasiewicz inequality*]{}, Math. Oper. Res., 35 (2010), pp. 438–457, <https://doi.org/10.1287/moor.1100.0449>.

, [*Dictionary learning for sparse coding: algorithms and convergence analysis*]{}, IEEE Trans. Pattern Anal. Mach. Intell., 38 (2016), pp. 1356–1369, <https://doi.org/10.1109/TPAMI.2015.2487966>.

, [*Convergence analysis for iterative data-driven tight frame construction scheme*]{}, Appl. Comput. Harmon. Anal., 38 (2015), pp. 510–523, <https://doi.org/10.1016/j.acha.2014.06.007>.

, [*A second-order model for image denoising*]{}, Set-Valued Var. Anal., 18 (2010), pp. 277–306, <https://doi.org/10.1007/s11228-010-0156-6>.

, [*Atomic norm denoising with applications to line spectral estimation*]{}, IEEE Trans. Signal Process., 61 (2013), pp. 5987–5999, <https://doi.org/10.1109/TSP.2013.2273443>.

, [ *Sparse sampling of signal innovations*]{}, IEEE Signal Process. Mag., 25 (2008), pp. 31–40, <https://doi.org/10.1109/MSP.2007.914998>.

, [*Proximal alternating linearized minimization for nonconvex and nonsmooth problems*]{}, Math. Program., 146 (2014), pp. 459–494, <https://doi.org/10.1007/s10107-013-0701-9>.

, [*Total generalized variation*]{}, SIAM J. Imaging Sci., 3 (2010), pp. 492–526, <https://doi.org/10.1137/090769521>.

, [*A framelet-based image inpainting algorithm*]{}, Appl. Comput. Harmon. Anal., 24 (2008), pp. 131–149, <https://doi.org/10.1016/j.acha.2007.10.002>.

, [*Image restoration: total variation, wavelet frames, and beyond*]{}, J. Amer. Math. Soc., 25 (2012), pp. 1033–1089, <https://doi.org/10.1090/S0894-0347-2012-00740-1>.

, [*Data-driven tight frame construction and image denoising*]{}, Appl. Comput. Harmon. Anal., 37 (2014), pp. 89–105, <https://doi.org/10.1016/j.acha.2013.10.001>.

, [*Split [B]{}regman methods and frame based image restoration*]{}, Multiscale Model. Simul., 8 (2009/10), pp. 337–369, <https://doi.org/10.1137/090753504>.

, [*Robust recovery of complex exponential signals from random [G]{}aussian projections via low rank [H]{}ankel matrix reconstruction*]{}, Appl. Comput. Harmon. Anal., 41 (2016), pp. 470–490, <https://doi.org/10.1016/j.acha.2016.02.003>.

, [*Spectral compressed sensing via projected gradient descent*]{}, SIAM J. Optim., 28 (2018), pp. 2625–2653, <https://doi.org/10.1137/17M1141394>.

, [*Fast and provable algorithms for spectrally sparse signal reconstruction via low-rank [H]{}ankel matrix completion*]{}, Appl. Comput. Harmon. Anal., 46 (2019), pp. 94–121, <https://doi.org/10.1016/j.acha.2017.04.004>.

, [*Super-resolution from noisy data*]{}, J. Fourier Anal. Appl., 19 (2013), pp. 1229–1254, <https://doi.org/10.1007/s00041-013-9292-3>.

, [*Towards a mathematical theory of super-resolution*]{}, Comm. Pure Appl. Math., 67 (2014), pp. 906–956, <https://doi.org/10.1002/cpa.21455>.

, [*Exact matrix completion via convex optimization*]{}, Found. Comput. Math., 9 (2009), pp. 717–772, <https://doi.org/10.1007/s10208-009-9045-5>.

, [*Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information*]{}, IEEE Trans. Inform. Theory, 52 (2006), pp. 489–509, <https://doi.org/10.1109/TIT.2005.862083>.

, [*Wavelet algorithms for high-resolution image reconstruction*]{}, SIAM J. Sci. Comput., 24 (2003), pp. 1408–1432, <https://doi.org/10.1137/S1064827500383123>.

, [*Robust spectral compressed sensing via structured matrix completion*]{}, IEEE Trans. Inform. Theory, 60 (2014), pp. 6576–6601, <https://doi.org/10.1109/TIT.2014.2343623>.

, [*Sensitivity to basis mismatch in compressed sensing*]{}, in 2010 IEEE International Conference on Acoustics, Speech and Signal Processing, March 2010, pp. 3930–3933, <https://doi.org/10.1109/ICASSP.2010.5495800>.

, [*P[ET]{}-[MRI]{} joint reconstruction by joint sparsity based tight frame regularization*]{}, SIAM J. Imaging Sci., 11 (2018), pp. 1179–1204, <https://doi.org/10.1137/17M1131453>.

, [*M[RA]{}-based wavelet frames and applications*]{}, in Mathematics in Image Processing, vol. 19 of IAS/Park City Math. Ser., Amer. Math. Soc., Providence, RI, 2013, pp. 9–158.

, [*Image restoration: a data-driven perspective*]{}, in Proceedings of the 8th [I]{}nternational [C]{}ongress on [I]{}ndustrial and [A]{}pplied [M]{}athematics, Higher Ed. Press, Beijing, 2015, pp. 65–108.

, [*Image restoration: a general wavelet frame based model and its asymptotic analysis*]{}, SIAM J. Math. Anal., 49 (2017), pp. 421–445, <https://doi.org/10.1137/16M1064969>.

, [*Compressed sensing*]{}, IEEE Trans. Inform. Theory., 52 (2006), pp. 1289–1306, <https://doi.org/10.1109/TIT.2006.871582>.

, [*Hankel matrix rank minimization with applications to system identification and realization*]{}, SIAM J. Matrix Anal. Appl., 34 (2013), pp. 946–977, <https://doi.org/10.1137/110853996>.

, [*Real Analysis: Modern Techniques and Their Applications*]{}, Pure and Appl. Math., John Wiley & Sons Inc., New York, 2nd ed., 1999.

, [*Low-rank matrix recovery via iteratively reweighted least squares minimization*]{}, SIAM J. Optim., 21 (2011), pp. 1614–1640, <https://doi.org/10.1137/100811404>.

, [*The split [B]{}regman method for [$L1$]{}-regularized problems*]{}, SIAM J. Imaging Sci., 2 (2009), pp. 323–343, <https://doi.org/10.1137/080725891>.

, [ *Realistic analytical phantoms for parallel magnetic resonance imaging*]{}, IEEE Trans. Med. Imag., 31 (2012), pp. 626–636, <https://doi.org/10.1109/TMI.2011.2174158>.

, [ *Magnetic Resonance Imaging : Physical Principles and Sequence Design*]{}, Wiley, 1st ed., June 1999.

, [*Adaptive multiresolution analysis structures and shearlet systems*]{}, SIAM J. Numer. Anal., 49 (2011), pp. 1921–1946, <https://doi.org/10.1137/090780912>.

, [*Deep learning for undersampled [MRI]{} reconstruction*]{}, Phys. Med. Biol., 63 (2018), p. 135007, <https://doi.org/10.1088/1361-6560/aac71a>.

, [*Modern Spectral Estimation: Theory and Application*]{}, Prentice-Hall Signal Processing Series: Advanced monographs, PTR Prentice Hall, 1988.

, [*Spectrum analysis?a modern perspective*]{}, Proceedings of the IEEE, 69 (1981), pp. 1380–1419, <https://doi.org/10.1109/PROC.1981.12184>.

, [*On gradients of functions definable in o-minimal structures*]{}, Ann. Inst. Fourier (Grenoble), 48 (1998), pp. 769–783, <http://www.numdam.org/item?id=AIF_1998__48_3_769_0>.

, [ *Balanced sparse model for tight frames in compressed sensing magnetic resonance imaging*]{}, PLOS ONE, 10 (2015), pp. 1–19, <https://doi.org/10.1371/journal.pone.0119584>.

, [*On semi-analytic and subanalytic geometry*]{}, in Panoramas of Mathematics ([W]{}arsaw, 1992/1994), vol. 34 of Banach Center Publ., Polish Acad. Sci. Inst. Math., Warsaw, 1995, pp. 89–104.

, [*Sparse [MRI]{}: the application of compressed sensing for rapid [MR]{} imaging*]{}, Magn. Reson. Med., 58 (2007), pp. 1182–1195, <https://doi.org/10.1002/mrm.21391>.

, [*Iterative reweighted algorithms for matrix rank minimization*]{}, J. Mach. Learn. Res., 13 (2012), pp. 3441–3473.

, [*Convex recovery of continuous domain piecewise constant images from nonuniform [F]{}ourier samples*]{}, IEEE Trans. Signal Process., 66 (2018), pp. 236–250, <https://doi.org/10.1109/TSP.2017.2750111>.

, [*Recovery of piecewise smooth images from few fourier samples*]{}, in 2015 International Conference on Sampling Theory and Applications (SampTA), May 2015, pp. 543–547, <https://doi.org/10.1109/SAMPTA.2015.7148950>.

, [*Super-resolution [MRI]{} using finite rate of innovation curves*]{}, in 2015 IEEE 12th International Symposium on Biomedical Imaging (ISBI), April 2015, pp. 1248–1251, <https://doi.org/10.1109/ISBI.2015.7164100>.

, [*Off-the-grid recovery of piecewise constant images from few [F]{}ourier samples*]{}, SIAM J. Imaging Sci., 9 (2016), pp. 1004–1041, <https://doi.org/10.1137/15M1042280>.

, [*A fast algorithm for convolutional structured low-rank matrix recovery*]{}, IEEE Trans. Comput. Imag., 3 (2017), pp. 535–550, <https://doi.org/10.1109/TCI.2017.2721819>.

, [*Sampling curves with finite rate of innovation*]{}, IEEE Trans. Signal Process., 62 (2014), pp. 458–471, <https://doi.org/10.1109/TSP.2013.2292033>.

, [*A combined first and second order variational approach for image reconstruction*]{}, J. Math. Imaging Vision, 48 (2014), pp. 308–338, <https://doi.org/10.1007/s10851-013-0445-4>.

, [*Essai expérimental et analytique: sur les lois de la dilatabilite des fluids elastiques et sur celles della force expansive de la vapeur de l’eau et de la vapeur de l’alkool, à differentes temperatures*]{}, J. Ec. Polytech. Paris, 1 (1795), pp. 24–76.

, [*Mr image reconstruction from highly undersampled k-space data by dictionary learning*]{}, IEEE Trans. Med. Imag., 30 (2011), pp. 1028–1041, <https://doi.org/10.1109/TMI.2010.2090538>.

, [ *Super-resolution in magnetic resonance imaging: [A]{} review*]{}, Concept. Magn. Reson. A, 40A (2012), pp. 306–325, <https://doi.org/10.1002/cmr.a.21249>.

, [*Affine systems in [$L_2(\bold R^d)$]{}: the analysis of the analysis operator*]{}, J. Funct. Anal., 148 (1997), pp. 408–447, <https://doi.org/10.1006/jfan.1996.3079>.

, [*Wavelet frames and image restorations*]{}, in Proceedings of the [I]{}nternational [C]{}ongress of [M]{}athematicians. [V]{}olume [IV]{}, Hindustan Book Agency, New Delhi, 2010, pp. 2834–2863.

, [*Introduction to Spectral Analysis*]{}, Prentice Hall, 1997.

, [*Near minimax line spectral estimation*]{}, IEEE Trans. Inform. Theory, 61 (2015), pp. 499–512, <https://doi.org/10.1109/TIT.2014.2368122>.

, [*Data-driven tight frame for multi-channel images and its application to joint color-cepth image reconstruction*]{}, J. Oper. Res. Soc. China, 3 (2015), pp. 99–115, <https://doi.org/10.1007/s40305-015-0074-2>.

, [*Precise semidefinite programming formulation of atomic norm minimization for recovering d-dimensional (d ≥ 2) off-the-grid frequencies*]{}, in 2014 Information Theory and Applications Workshop (ITA), Feb 2014, pp. 1–4, <https://doi.org/10.1109/ITA.2014.6804267>.

, [*A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion*]{}, SIAM J. Imaging Sci., 6 (2013), pp. 1758–1789, <https://doi.org/10.1137/120887795>.

, [*Hankel matrix nuclear norm regularized tensor completion for [$N$]{}-dimensional exponential signals*]{}, IEEE Trans. Signal Process., 65 (2017), pp. 3702–3717, <https://doi.org/10.1109/TSP.2017.2695566>.

, [*C[T]{} image reconstruction by spatial-[R]{}adon domain data-driven tight frame regularization*]{}, SIAM J. Imaging Sci., 9 (2016), pp. 1063–1083, <https://doi.org/10.1137/16M105928X>.

, [*Coil compression for accelerated imaging with cartesian sampling*]{}, Magn. Reson. Med., 69 (2013), pp. 571–582, <https://doi.org/10.1002/mrm.24267>, <https://arxiv.org/abs/https://onlinelibrary.wiley.com/doi/pdf/10.1002/mrm.24267>.

[^1]: Department of Mathematics, Hong Kong University of Science and Technology, Clearwater Bay, Kowloon, Hong Kong, China ().

[^2]: Corresponding Author. School of Mathematical Sciences, Tongji University, 1239 Siping Road, Shanghai, 200092 China ().

[^3]: Corresponding Author. School of Data Science, Fudan University, 220 Handan Road, Shanghai, 200433 China ().

[^4]: Submitted to the editors DATE.
