text:
\section{Introduction}\label{Introduction}

Magnetic resonance imaging (MRI) is one of the most widely used medical imaging modality in clinical diagnosis \cite{E.M.Haacke1999}. It is non-radioactive, non-invasive, and has excellent soft tissue contrasts such as T1 and T2 with high spatial resolution \cite{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 \cite{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 (motion artifacts due to the patient's respiratory motion which can generate the outlier of the k-space data) constraints \cite{C.M.Hyun2018,K.H.Jin2017,Y.Liu2015,M.Lustig2007}, there has been increasing demand for methods which can reduce the amount of acquired data without degrading the image quality \cite{M.Lustig2007}.

When the k-space data is undersampled, the Nyquist sampling criterion is violated, and the \emph{direct restoration} using this undersampled k-space data will lead to the aliasing in the reconstructed image \cite{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 \cite{E.J.Candes2006,Donoho2006,C.M.Hyun2018,M.Lustig2007}. More precisely, a typical $\ell_1$ norm based ``on-the-grid'' CS-MRI model restores the MR image on the grid, by regularizing smooth image components while preserving image singularities such as edges, ridges, and corners.

Even though the conventional on-the-grid approaches have shown strong ability to reduce the data acquisition time and thus received a lot of attention over the past few years \cite{Y.Liu2015}, they have 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 \emph{continuous domain}, under the assumption of a single-coil MRI \cite{E.M.Haacke1999}. This means that, the sparse regularization based image restoration will work well when the singularities of the image are well aligned with the grid. However, even in the case of piecewise constant function whose (distributional) gradient is sparse in the continuous domain, its singularities may not necessarily agree with the discrete grid, leading to the problem of \emph{basis mismatch} \cite{Y.Chi2010,G.Ongie2016,J.Ying2017}. Such a basis mismatch between the true singularities and the discrete grid may result in the loss of sparse structure of the image, and thus degrades the restoration quality \cite{J.Ying2017}.

\subsection{Continuous domain regularization for image restoration}

The continuous domain regularization is emerging as a powerful alternative to the discrete domain sparse regularization \cite{B.N.Bhaskar2013,E.J.Candes2014,Y.Chen2014,G.Ongie2018}. By such an ``off-the-grid'' scheme, the regularization on the discrete data exploits the sparsity in the continuous domain to alleviate the basis mismatch, which is especially attractive in signal restorations from partial measurements \cite{G.Ongie2018}. One of the most successful examples is the finite-rate-of-innovation (FRI) framework which corresponds a superposition of a few sinusoids to a low rank Hankel matrix \cite{T.Blu2008}. Similarly, we can adopt the so-called \emph{structured low rank matrix completion (SLMC)} \cite{Y.Chen2014,G.Ongie2017,J.C.Ye2017} for the CS restoration based on the continuous domain regularization. The SLMC approach can also be used to restore Fourier samples of a one dimensional piecewise constant signal \cite{T.Blu2008,M.Vetterli2002} by considering the Fourier samples of a derivative. However, while this framework works well in the case of isolated singularities \cite{E.J.Candes2013,W.Xu2014}, the extension to the (piecewise smooth) image restoration is not straightforward \cite{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 makes it challenging to construct a structured low rank matrix from the Fourier samples.

Recently, there are several extensions of the FRI framework to two dimensional continuous domain regularization for image restoration. Such extensions include the piecewise holomorphic complex image restoration \cite{H.Pan2014} and the piecewise constant real image restoration \cite{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 \emph{annihilating polynomial}. Indeed, given that the edges lie in the zero level set of an annihilating polynomial, the structured matrix (Hankel/Toeplitz matrix) constructed from the Fourier samples of piecewise constant images becomes low rank \cite{G.Ongie2018,G.Ongie2015,G.Ongie2016}. Hence, in CS-MRI, we can first restore the fully sampled k-space data (the discrete sample of the Fourier transform of a piecewise constant function) via the SLMC, and then restore MR image by the inverse DFT \cite{Haldar2014,K.H.Jin2016}.

Though the rank minimization problem is an NP-hard problem in general \cite{E.J.Candes2009}, numerous tractable relaxation approaches have been proposed for the SLMC. One of them is the convex nuclear norm relaxation method (e.g. \cite{J.F.Cai2016,M.Fazel2013}), together with the theoretical restoration guarantees (e.g. \cite{G.Ongie2018}). In addition, the iterative reweighted least squares (IRLS) approaches for the Schatten $p$-norm minimization are proposed in \cite{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 \cite{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.

\subsection{Motivations and contributions of our approach}

Our proposed off-the-grid approach for the CS-MRI restoration is based on the fact that the SVD of a Hankel matrix induces tight frame filter banks due to its underlying convolutional structure. This means that, if we can associate a signal with a low rank Hankel matrix, its right singular vectors form tight frame filter banks which allow us to represent the signal with sparse canonical coefficients. Hence, instead of directly minimizing the rank of a two-fold Hankel matrix, we can develop the sparse regularization model via data driven tight frames \cite{J.F.Cai2014} for the CS-MRI restoration. More specifically, using the $\ell_0$ norm of the frame coefficients, we propose the \emph{balanced approach} \cite{J.F.Cai2008,R.H.Chan2003} as a relaxation of the low rank two-fold Hankel matrix approximation which is based on the framework of annihilating filters. Finally, the numerical experiments show that our data driven tight frame approach outperforms the structured low rank matrix approaches (and their relaxations) as well as the existing on-the-grid approaches, leading to the state-of-the-art performance.

Note that the data driven tight frame as well as the dictionary learning has been discussed for the various image restoration tasks. More precisely, the data driven tight frame is first used for the image denoising \cite{J.F.Cai2014}, the multi-channel image denoising and inpainting \cite{J.Wang2015}, the joint spatial-Radon domain regularization for the sparse view CT \cite{R.Zhan2016}, and the PET-MRI joint reconstruction \cite{J.K.Choi2018}. In addition, the authors in \cite{S.Ravishankar2015} explored the general framework of the dictionary learning based image restoration model with the convergence analysis. However, while these aforementioned approaches mainly concern the sparse approximation of the discrete images, to the best of our knowledge, our approach for the first time applies the data driven tight frame regularization to the MRI restoration based on a continuous image model.

Finally, we would like to mention that the proposed off-the-grid approach for the CS-MRI restoration is not limited to the two dimensional single-coil MRI restoration. Even though further numerical studies are required, it is not hard to generalize the proposed model to the three dimensional cases. In addition, since the multiplication in the image domain becomes the convolution in the frequency domain via the Fourier transform, we can extend our approach to a SENSE setup \cite{K.P.Pruessmann1999}, a multiple-coil MRI restoration whose inverse problem corresponds to the multi-channel deconvolution in the frequency domain with a known coil sensitivity map.

\subsection{Organization and notation of paper}

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

\begin{figure}[tp!]
\centering
\includegraphics[width=0.95\textwidth]{HankelIllustrate.pdf}\vspace{-0.20cm}
\caption{Illustration of generating a Hankel matrix ${\boldsymbol \mH}\boldsymbol{v}$ from $\boldsymbol{v}\in{\mathcal V}_2$ and a two-fold Hankel matrix ${\boldsymbol \mH}\mathbf{w}$ from $\mathbf{w}=\left(\mathbf{w}_1,\mathbf{w}_2\right)\in{\mathcal V}_2\times{\mathcal V}_2$. Here, $K_1=K_2=2$.}\label{BlockHankelIllustrate}
\end{figure}

Throughout this paper, all two dimensional discrete images will be denoted by the bold faced lower case letters while the functions in a two dimensional continuous domain will be denoted by regular characters. Note that a two dimensional image can also be identified with a vector whenever convenient. All matrices will be denoted by the bold faced upper case letters, and the $j$th row and the $k$th column of a matrix ${\boldsymbol Z}$ will be denoted by ${\boldsymbol Z}^{(j,:)}$ and ${\boldsymbol Z}^{(:,k)}$, respectively. We denote by
\begin{align}\label{ImageGrid}
{\mathbb O}=\left\{-\lfloor N/2\rfloor,\cdots,\lfloor N/2-1\rfloor\right\}^2,
\end{align}
with $N\in{\mathbb N}$, the set of $N\times N$ grid, and the space of complex valued functions on ${\mathbb O}$ is denoted by ${\mathcal V}_2\simeq{\Bbb C}^{|{\mathbb O}|}$. Given two rectangular grids ${\mathbb K}$ and ${\mathbb M}$, the contraction ${\mathbb K}:{\mathbb M}$ is defined as
\begin{align*}
{\mathbb K}:{\mathbb M}=\left\{{\boldsymbol k}\in{\mathbb K}:{\boldsymbol k}-{\boldsymbol j}\in{\mathbb K}~\text{for all}~{\boldsymbol j}\in{\mathbb M}\right\}.
\end{align*}
Finally, operators are denoted as the bold faced calligraphic letters. For instance, let $\boldsymbol{v}\in{\mathcal V}_2$ and let ${\mathbb K}$ be a rectangular $K_1\times K_2$ grid. The corresponding Hankel matrix ${\boldsymbol \mH}\boldsymbol{v}$ is an $M_1\times M_2$ matrix ($M_1=|{\mathbb O}:{\mathbb K}|$ and $M_2=|{\mathbb K}|$) generated by concatenating $K_1\times K_2$ patches of $\boldsymbol{v}$ into row vectors. In the sense of two dimensional multi-indices, we can formally write ${\boldsymbol \mH}\boldsymbol{v}$ as
\begin{align*}
\left({\boldsymbol \mH}\boldsymbol{v}\right)[{\boldsymbol k},{\boldsymbol l}]=\boldsymbol{v}[{\boldsymbol k}+{\boldsymbol l}],~~~~~~{\boldsymbol k}\in{\mathbb O}:{\mathbb K}~~~\text{and}~~~{\boldsymbol l}\in{\mathbb K}.
\end{align*}
With a slight abuse of notation, we also use
\begin{align}\label{TwoFoldHankel}
{\boldsymbol \mH}\mathbf{w}=\left[\begin{array}{cc}
\left({\boldsymbol \mH}\mathbf{w}_1\right)^T&\left({\boldsymbol \mH}\mathbf{w}_2\right)^T
\end{array}\right]^T,
\end{align}
to denote the $2M_1\times M_2$ two-fold Hankel matrix constructed from $\mathbf{w}=\left(\mathbf{w}_1,\mathbf{w}_2\right)\in{\mathcal V}_2\times{\mathcal V}_2$; see \cref{BlockHankelIllustrate} for an illustration.

\section{Preliminaries on tight wavelet frames}\label{TightFrames}

Here we provide a brief introduction on tight wavelet frames and data driven tight frames. Interested readers may consult \cite{J.F.Cai2008,B.Dong2013,B.Dong2015,A.Ron1997,Shen2010} for detailed surveys on tight wavelet frames, and \cite{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 ${\mathscr H}$ a Hilbert space and let $\langle\cdot,\cdot\rangle$ be the inner product defined on ${\mathscr H}$. A countable set $\{\boldsymbol{\varphi}_n:n\in{\mathbb Z}\}\subseteq{\mathscr H}$ is called a tight frame on ${\mathscr H}$ if
\begin{align}\label{TightFrame}
\|\boldsymbol{u}\|^2=\sum_{n\in{\mathbb Z}}|\langle\boldsymbol{u},\boldsymbol{\varphi}_n\rangle|^2~~~~~\text{for all}~~~\boldsymbol{u}\in{\mathscr H}.
\end{align}
Given $\{\boldsymbol{\varphi}_n:n\in{\mathbb Z}\}\subseteq{\mathscr H}$, we define the analysis operator ${\boldsymbol \mW}:{\mathscr H}\to\ell_2({\mathbb Z})$ as
\begin{align*}
\boldsymbol{u}\in{\mathscr H}\mapsto {\boldsymbol \mW}\boldsymbol{u}=\{\langle\boldsymbol{u},\boldsymbol{\varphi}_n\rangle:n\in{\mathbb Z}\}\in\ell_2({\mathbb Z}).
\end{align*}
The synthesis operator ${\boldsymbol \mW}^T:\ell_2({\mathbb Z})\to{\mathscr H}$ is defined as the adjoint of ${\boldsymbol \mW}$:
\begin{align*}
{\boldsymbol c}\in\ell_2({\mathbb Z})\mapsto{\boldsymbol \mW}^T{\boldsymbol c}=\sum_{n\in{\mathbb Z}}{\boldsymbol c}[n]\boldsymbol{\varphi}_n\in{\mathscr H}.
\end{align*}
Then $\{\boldsymbol{\varphi}_n:n\in{\mathbb Z}\}$ is a tight frame on ${\mathscr H}$ if and only if ${\boldsymbol \mW}^T{\boldsymbol \mW}={\boldsymbol \mI}$ where ${\boldsymbol \mI}$ is the identity on ${\mathscr H}$. It follows that, for a given tight frame $\{\boldsymbol{\varphi}_n:n\in{\mathbb Z}\}$, we have the following canonical expression:
\begin{align*}
\boldsymbol{u}=\sum_{n\in{\mathbb Z}}\langle\boldsymbol{u},\boldsymbol{\varphi}_n\rangle\boldsymbol{\varphi}_n,
\end{align*}
with ${\boldsymbol \mW}\boldsymbol{u}=\{\langle\boldsymbol{u},\boldsymbol{\varphi}_n\rangle:n\in{\mathbb 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 $\|\boldsymbol{\varphi}_n\|=1$ for all $n\in{\mathbb Z}$.

One of the most widely used class of tight frames is the discrete wavelet frame generated by a set of finitely supported filters $\{\mathbf{a}_1,\cdots,\mathbf{a}_m\}$. Throughout this paper, we only discuss the two dimensional undecimated wavelet frames on $\ell_2({\mathbb Z}^2)$, but note that it is not difficult to extend to $\ell_2({\mathbb Z}^d)$ with $d\geq3$. For $\mathbf{a}\in\ell_1({\mathbb Z}^2)$, define a convolution operator ${\boldsymbol \mS}_{\mathbf{a}}:\ell_2({\mathbb Z}^2)\to\ell_2({\mathbb Z}^2)$ by
\begin{align}\label{DiscreteConv}
({\boldsymbol \mS}_{\mathbf{a}}\boldsymbol{u})[{\boldsymbol k}]=(\mathbf{a}\ast\boldsymbol{u})[{\boldsymbol k}]=\sum_{{\boldsymbol l}\in{\mathbb Z}^2}\mathbf{a}[{\boldsymbol k}-{\boldsymbol l}]\boldsymbol{u}[{\boldsymbol l}]~~~~~\text{for}~~~\boldsymbol{u}\in\ell_2({\mathbb Z}^2).
\end{align}
Given a set of finitely supported filters $\{\mathbf{a}_1,\cdots,\mathbf{a}_m\}$, define the analysis operator ${\boldsymbol \mW}$ and the synthesis operator ${\boldsymbol \mW}^T$ by
\begin{align}
{\boldsymbol \mW}&=\left[{\boldsymbol \mS}_{\mathbf{a}_1[-\cdot]}^T,{\boldsymbol \mS}_{\mathbf{a}_2[-\cdot]}^T,\cdots,{\boldsymbol \mS}_{\mathbf{a}_m[-\cdot]}^T\right]^T,\label{AnalConv}\\
{\boldsymbol \mW}^T&=\big[{\boldsymbol \mS}_{\mathbf{a}_1},{\boldsymbol \mS}_{\mathbf{a}_2},\cdots,{\boldsymbol \mS}_{\mathbf{a}_m}\big],\label{SyntConv}
\end{align}
respectively. Then, the direct computation can show that the rows of ${\boldsymbol \mW}$ form a tight frame on $\ell_2({\mathbb Z}^2)$ (i.e. ${\boldsymbol \mW}^T{\boldsymbol \mW}={\boldsymbol \mI}$) if and only if the filters $\{\mathbf{a}_1,\cdots,\mathbf{a}_m\}$ satisfy
\begin{align}\label{SomeUEP}
\sum_{j=1}^m\sum_{{\boldsymbol l}\in{\mathbb Z}^2}\mathbf{a}_j[{\boldsymbol k}+{\boldsymbol l}]\mathbf{a}_j[{\boldsymbol l}]=\boldsymbol{\delta}[{\boldsymbol k}]=\left\{\begin{array}{cl}
1~&\text{if}~{\boldsymbol k}={\mathbf 0},\\
0~&\text{if}~{\boldsymbol k}\neq{\mathbf 0},
\end{array}\right.
\end{align}
called the {\emph{unitary extension principle}} (UEP) condition \cite{B.Han2011}. Finally, for a two dimensional discrete image on the finite grid, ${\boldsymbol \mS}_{\mathbf{a}}$ in \cref{DiscreteConv} with a finitely supported $\mathbf{a}$ denotes the discrete convolution under the periodic boundary condition throughout this paper.

In the literature, wavelet frames are widely used for the sparse approximation of an image. To explore a better sparse approximation, the authors in \cite{J.F.Cai2014} proposed a \emph{data driven tight frame} approach, inspired by \cite{M.Aharon2006,S.Ravishankar2013}. Specifically, given a two dimensional image $\boldsymbol{u}\in{\mathcal V}_2$, a tight frame system ${\boldsymbol \mW}$ defined as in \cref{AnalConv}, which is generated by a set of finitely supported $p\times p$ filters $\{\mathbf{a}_1,\cdots,\mathbf{a}_{p^2}\}$ supported on a $p\times p$ grid ${\mathbb K}$ and satisfying \cref{SomeUEP}, is constructed via the following minimization
\begin{align}\label{DDTFModel}
\min_{{\boldsymbol c},{\boldsymbol \mW}}~\|{\boldsymbol c}-{\boldsymbol \mW}\boldsymbol{u}\|_2^2+\gamma^2\|{\boldsymbol c}\|_0~~~~~\text{subject to}~~~~{\boldsymbol \mW}^T{\boldsymbol \mW}={\boldsymbol \mI},
\end{align}
with the $\ell_0$ norm $\|{\boldsymbol c}\|_0$ encoding the number of nonzero entries in the coefficient vector ${\boldsymbol c}$. Note that in wavelet domain, there is usually a set of low pass coefficients that are not sparse. In fact, through the relaxation of \cref{DDTFModel} which will be discussed below, we can learn a tight frame system ${\boldsymbol \mW}$ by fixing a low pass filter and penalizing the $\ell_0$ norm of high pass coefficients only, following \cite{C.Bao2013}. However, since the main focus is to establish a relation between the tight frame filter banks and the Hankel matrix, we forgo further discussions on such a detailed construction.

To solve \cref{DDTFModel}, the authors in \cite{J.F.Cai2014} presented that if ${\boldsymbol A}\in{\Bbb R}^{p^2\times p^2}$ satisfies ${\boldsymbol A}\bsA^T=p^{-2}{\boldsymbol I}$, after reformulating into $p\times p$ filters supported on ${\mathbb K}$, its column vectors generate filters satisfying the UEP condition \cref{SomeUEP}. Hence, under such an assumption on the filters $\mathbf{a}_1,\cdots,\mathbf{a}_{p^2}$, we can rewrite \cref{DDTFModel} as
\begin{align}\label{DDTFModelModi}
\min_{{\boldsymbol C},{\boldsymbol A}}~\|{\boldsymbol C}-\left({\boldsymbol \mH}\boldsymbol{u}\right){\boldsymbol A}\|_F^2+\gamma^2\|{\boldsymbol C}\|_0~~~~~\text{subject to}~~~~{\boldsymbol A}\bsA^T=p^{-2}{\boldsymbol I},
\end{align}
to construct filter banks. Here, ${\boldsymbol \mH}\boldsymbol{u}\in{\Bbb R}^{(N-p+1)^2\times p^2}$ is a Hankel matrix generated by the $p\times p$ patches of $\boldsymbol{u}$, ${\boldsymbol C}$ is the corresponding matrix formulation of ${\boldsymbol c}$, and $\left\|\cdot\right\|_F$ is the Frobenius norm of a matrix. To solve \cref{DDTFModelModi}, the alternating minimization method with closed form solutions for each stage is presented in \cite{J.F.Cai2014}. In addition, the proximal alternating minimization (PAM) scheme with global convergence to critical points is proposed in \cite{C.Bao2015}.

\section{Proposed off-the-grid regularization model}\label{ModelandAlgorithm}

\subsection{Structured low rank matrix approaches in CS-MRI}\label{SRLACSMRI}

We begin with the brief review on the structured low rank matrix approaches in CS-MRI. Interested readers can refer to \cite{G.Ongie2016} and references therein for more detailed surveys. Throughout this paper, we only consider the two dimensional single-coil MRI reconstruction for the sake of simplicity, so that the k-space data ${\boldsymbol f}$ can be modeled as
\begin{align}\label{MRForward}
{\boldsymbol f}[{\boldsymbol k}]=\widehat{u}(L^{-1}{\boldsymbol k})=\int_{{\Bbb R}^2}u(\mathbf{x})e^{-2\pi i\mathbf{x}\cdot{\boldsymbol k}/L}\mathrm{d}\mathbf{x},~~~~~{\boldsymbol k}\in{\mathbb M},
\end{align}
where ${\mathbb M}$ is the sampling region, $L>0$ is the length of field-of-view (FOV), and $u\in L_1({\Bbb R}^2)$ is the proton spin density distribution in ${\Bbb R}^2$ \cite{E.M.Haacke1999}. We further assume that $u$ is supported on $[-L/2,L/2)^2$. This means that when ${\boldsymbol f}[{\boldsymbol k}]$ in \cref{MRForward} is available for ${\boldsymbol k}\in{\mathbb Z}^2$, we have
\begin{align*}
u(\mathbf{x})=L^{-2}\sum_{{\boldsymbol k}\in{\mathbb Z}^2}{\boldsymbol f}[{\boldsymbol k}]e^{2\pi i{\boldsymbol k}\cdot\mathbf{x}/L},~~~~~\mathbf{x}\in[-L/2,L/2)^2.
\end{align*}
In general, the k-space data can be sampled at ``arbitrary'' positions, i.e. ${\mathbb M}\subseteq{\Bbb R}^2$. However, to simplify the discussion, we assume that ${\mathbb M}\subseteq{\mathbb O}\subseteq{\mathbb Z}^2$ where ${\mathbb O}$ is defined as \cref{ImageGrid}, i.e. the discrete sampling on the Nyquist rate grid \cite{E.M.Haacke1999,T.H.Kim2019}.

The conventional CS-MRI approach aims to directly restore the MR image $\boldsymbol{u}$ on ${\mathbb O}$ from the k-space ${\boldsymbol f}$ undersampled on ${\mathbb M}$. Note that when ${\mathbb M}={\mathbb O}$, i.e. the fully sampled case, $\boldsymbol{u}$ is obtained from the inverse DFT
\begin{align}\label{MRInverse}
\boldsymbol{u}[\boldsymbol{p}]=\left({\boldsymbol \msF}^{-1}{\boldsymbol f}\right)[\boldsymbol{p}]=N^{-2}\sum_{{\boldsymbol k}\in{\mathbb O}}{\boldsymbol f}[{\boldsymbol k}]e^{2\pi i\boldsymbol{p}\cdot{\boldsymbol k}/N},~~~~~~~\boldsymbol{p}\in{\mathbb O}.
\end{align}
Hence, the conventional sparse regularization based CS-MRI approach
\begin{align}\label{CSApproach}
\min_{\boldsymbol{u}}~\frac{1}{2}\left\|{\boldsymbol \mR}_{{\mathbb M}}{\boldsymbol \msF}\boldsymbol{u}-{\boldsymbol f}\right\|_2^2+\gamma\left\|\boldsymbol{\Phi}\boldsymbol{u}\right\|_1
\end{align}
is based on the DFT
\begin{align}\label{DFT}
{\boldsymbol f}[{\boldsymbol k}]=\left({\boldsymbol \msF}\boldsymbol{u}\right)[{\boldsymbol k}]=\sum_{\boldsymbol{p}\in{\mathbb O}}\boldsymbol{u}[\boldsymbol{p}]e^{-2\pi i\boldsymbol{p}\cdot{\boldsymbol k}/N},~~~~~~~{\boldsymbol k}\in{\mathbb M}.
 \end{align}
Here, ${\boldsymbol \mR}_{{\mathbb M}}$ is a restriction onto ${\mathbb M}$, and $\boldsymbol{\Phi}$ is a sparsifying transformation. Note that \cref{CSApproach} is the on-the-grid scheme, as we promote the sparsity of the discrete image $\boldsymbol{u}$ on the grid ${\mathbb O}$. In other words, we can apply the sparse regularization effectively provided that the singularities of $u$ are well aligned with ${\mathbb O}$. However, since the k-space data actually comes from the Fourier transform of a continuous domain function $u$, together with \cref{MRForward} and ${\mathbb M}={\mathbb O}$, \cref{MRInverse} can be rewritten as
\begin{align*}
\boldsymbol{u}[\boldsymbol{p}]=\left(N^{-1}L\right)^2\left(D_N\ast u\right)(N^{-1}L\boldsymbol{p}),~~~~~\boldsymbol{p}\in{\mathbb O},
\end{align*}
where $D_N(\mathbf{x})$ is the Dirichlet kernel defined as
\begin{align*}
D_N(\mathbf{x})=L^{-2}\sum_{{\boldsymbol k}\in{\mathbb O}}e^{2\pi i{\boldsymbol k}\cdot\mathbf{x}/L}.
\end{align*}
In other words, since the MR image is proportional to a discrete sampling of $D_N\ast u$ on $N^{-1}L{\mathbb O}$ rather than that of $u$ \cite{E.M.Haacke1999}, there exists a basis mismatch between the singularities of $u$ in the continuous domain and the discrete grid ${\mathbb O}$. 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 $\boldsymbol{v}$ from its undersampled version ${\boldsymbol f}$ modeled as in \cref{MRForward}. Formally, we can write the problem as follows:
\begin{align}\label{OurProblem}
\text{find}~~\boldsymbol{v}={\mathscr F}(u)\big|_{L^{-1}{\mathbb O}}~~\text{subject to}~~{\boldsymbol \mR}_{{\mathbb M}}\boldsymbol{v}={\boldsymbol f},
\end{align}
and then compute $\boldsymbol{u}={\boldsymbol \msF}^{-1}\boldsymbol{v}$ by using the restored k-space data $\boldsymbol{v}$.

\begin{figure}[tp!]
\centering
\includegraphics[width=0.95\textwidth]{CSMRIComparisons.pdf}\vspace{-0.20cm}
\caption{Comparison between the conventional CS-MRI and the off-the-grid CS-MRI approach.}\label{CSMRIComparisons}
\end{figure}

Generally, it is impossible to directly solve \cref{OurProblem} without any further information. Nevertheless, by imposing the prior information of $u$ (in the continuous domain) into \cref{OurProblem}, we can exploit the continuous domain regularization, and reduce the basis mismatch; see \cref{CSMRIComparisons} for the schematic comparison. For this purpose, we follow \cite{G.Ongie2016} and consider the following piecewise constant function model of the proton density $u$:
\begin{align}\label{uModel}
u(\mathbf{x})=\sum_{j=1}^J\alpha_j1_{\Omega_j}(\mathbf{x}),~~~~~~~\mathbf{x}\in{\Bbb R}^2,
\end{align}
where $\alpha_j\in{\Bbb C}$ and $1_{\Omega}$ denotes the characteristic function on a set $\Omega$: $1_{\Omega}(\mathbf{x})=1$ if $\mathbf{x}\in\Omega$, and $0$ otherwise. We further assume that each $\Omega_j$ lies in $[-L/2,L/2)^2$, and \cref{uModel} is expressed with the smallest number of characteristic functions such that $\Omega_j$'s are pairwise disjoint. Then the discontinuities of $u$ agrees with $\Gamma=\bigcup_{j=1}^J\partial\Omega_j$, called the \emph{edge set} of $u$ \cite{G.Ongie2016}. Under this setting, the authors in \cite{G.Ongie2015,G.Ongie2016} derived the following \emph{linear annihilation relation}.

\begin{proposition}\label{Prop1} Let $u(\mathbf{x})$ be a piecewise constant function defined as \cref{uModel}. Assume that there exists a finite (rectangular and symmetric) index set ${\mathbb K}$ such that
\begin{align}\label{MajorAssumption}
\Gamma\subseteq\left\{\mathbf{x}\in{\Bbb R}^2:\varphi(\mathbf{x})=0\right\}~~~~~\text{with}~~\varphi(\mathbf{x})=\sum_{{\boldsymbol k}\in{\mathbb K}}\mathbf{a}[{\boldsymbol k}]e^{2\pi i{\boldsymbol k}\cdot\mathbf{x}/L}.
\end{align}
Then for such $\varphi$, the Fourier transform of $\nabla u=\left(\partial_1u,\partial_2u\right)$ is annihilated by the convolution with the Fourier transform of $\varphi$. More precisely, we have
\begin{align}\label{FouGradConvSemiDiscrete}
\left({\mathscr F}(\nabla u)\ast\widehat{\varphi}\right)({\boldsymbol{\xi}})=\sum_{{\boldsymbol k}\in{\mathbb K}}{\mathscr F}(\nabla u)({\boldsymbol{\xi}}-L^{-1}{\boldsymbol k})\mathbf{a}[{\boldsymbol k}]=0,~~~~~~{\boldsymbol{\xi}}\in{\Bbb R}^2,
\end{align}
where the convolution acts on ${\mathscr F}(\partial_1 u)$ and ${\mathscr F}(\partial_2 u)$ separately.
\end{proposition}

The trigonometric polynomial $\varphi$ in \cref{MajorAssumption} is called the \emph{annihilating polynomial}, $\mathbf{a}$ is called the \emph{annihilating filter} with support ${\mathbb K}$, and $\Gamma$ satisfying \cref{MajorAssumption} is called the \emph{trigonometric curves} \cite{G.Ongie2016}. It is proved in \cite{G.Ongie2016} that trigonometric curves can always be described as the zero set of a \emph{real valued} trigonometric polynomial, using the tool of algebraic geometries.

In MRI, the Fourier transform of $u$ is sampled on the grid $L^{-1}{\mathbb O}$, so \cref{FouGradConvSemiDiscrete} becomes the finite system of linear equations
\begin{align}\label{FouGradConvDiscrete}
\sum_{{\boldsymbol k}\in{\mathbb K}}{\mathscr F}(\nabla u)(L^{-1}({\boldsymbol l}+{\boldsymbol k}))\mathbf{a}[-{\boldsymbol k}]=0,~~~~~~~{\boldsymbol l}\in{\mathbb O}:{\mathbb K}.
\end{align}
In the matrix-vector multiplication form, we have
\begin{align}\label{RankDeficientHankel}
{\boldsymbol \mH}\left({\mathscr F}(\nabla u)\big|_{L^{-1}{\mathbb O}}\right)\mathbf{a}[-\cdot]={\mathbf 0}.
\end{align}
Moreover, if $\varphi$ in \cref{MajorAssumption} is the minimal polynomial with the support ${\mathbb K}$ (i.e. ${\mathbb K}$ has the smallest linear dimension) and ${\mathbb K}'$ is the assumed filter support strictly containing ${\mathbb K}$, then we have
\begin{align}\label{RankUB}
\mathrm{rank}\left({\boldsymbol \mH}\left({\mathscr F}(\nabla u)\big|_{L^{-1}{\mathbb O}}\right)\right)\leq|{\mathbb K}'|-|{\mathbb K}':{\mathbb K}|.
\end{align}
Roughly speaking, \cref{RankUB} means that if $\varphi$ in \cref{MajorAssumption} is the minimal polynomial with $\mathbf{a}$ satisfying \cref{FouGradConvDiscrete,RankDeficientHankel}, so does the translation $\mathbf{a}[\cdot-\boldsymbol{m}]$ for all $\boldsymbol{m}\in{\mathbb K}'\setminus\left({\mathbb K}':{\mathbb K}\right)$, or equivalently, $e^{2\pi i\boldsymbol{m}\cdot\mathbf{x}}\varphi(\mathbf{x})$ is also an annihilating polynomial. See \cite{G.Ongie2015,G.Ongie2016} for details.

Note that \cref{RankUB} means the two-fold Hankel matrix constructed by ${\mathscr F}(\nabla u)\big|_{L^{-1}{\mathbb O}}$ is \emph{rank deficient} whenever the grid ${\mathbb O}$ in \cref{ImageGrid} is large enough to satisfy the conditions in \cite[Proposition 5.1]{G.Ongie2016}. Based on this observation, we can formulate \cref{OurProblem} as the following low rank two-fold Hankel matrix completion \cite{Y.Chen2014,G.Ongie2017,J.C.Ye2017}
\begin{align}\label{LRSMCompletionRevist}
\min~~\mathrm{rank}\left({\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)\right)~~~~~~~\text{subject to}~~{\boldsymbol \mR}_{{\mathbb M}}\boldsymbol{v}={\boldsymbol f}.
\end{align}
Here, ${\boldsymbol\Lambda}\boldsymbol{v}=\left({\boldsymbol\Lambda}_1\boldsymbol{v},{\boldsymbol\Lambda}_2\boldsymbol{v}\right)\in{\mathcal V}_2\times{\mathcal V}_2$ and ${\boldsymbol\Lambda}\in{\Bbb C}^{2|{\mathbb O}|\times|{\mathbb O}|}$ is a weight matrix defined as
\begin{align}\label{Lambda}
{\boldsymbol\Lambda}=\left[\begin{array}{cc}
{\boldsymbol\Lambda}_1&{\boldsymbol\Lambda}_2
\end{array}\right]^T=\left[\begin{array}{cc}
\diag\left(2\pi ik_1/L\right)&\diag\left(2\pi ik_2/L\right)
\end{array}\right]_{{\boldsymbol k}=(k_1,k_2)\in{\mathbb O}}^T,
\end{align}
which is derived from the sample of ${\mathscr F}\left(\nabla u\right)({\boldsymbol{\xi}})=2\pi i{\boldsymbol{\xi}}\widehat{u}({\boldsymbol{\xi}})=\left(2\pi i\xi_1\widehat{u}({\boldsymbol{\xi}}),2\pi i\xi_2\widehat{u}({\boldsymbol{\xi}})\right)$ on $L^{-1}{\mathbb O}$. Hence, the piecewise constant property of $u$ in the continuous domain can be transformed to the low rank property of ${\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)$.

\subsection{Proposed CS-MRI model}\label{ProposedModelandAlgorithm}

Denote by $\boldsymbol{v}\in{\mathcal V}_2$ the k-space data on the grid ${\mathbb O}$, which is to be restored from a given data ${\boldsymbol f}$ undersampled on ${\mathbb M}$. Consider a two-fold Hankel matrix ${\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)\in{\Bbb C}^{2M_1\times M_2}$ defined as \cref{TwoFoldHankel} with ${\boldsymbol\Lambda}$ defined as \cref{Lambda}. To reconstruct $\boldsymbol{v}={\mathscr F}(u)\big|_{L^{-1}{\mathbb O}}$ corresponding to a piecewise constant function $u$ in \cref{uModel}, we assume that $\mathrm{rank}\left({\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)\right)=r\ll2M_1\wedge M_2$, following \cref{SRLACSMRI}. Considering its SVD, we present \cref{Th1}, which is the key idea of the proposed CS-MRI model. The proof can be found in \cref{ProofTh1}.

To do this, we let $\mathbf{a}_1,\cdots,\mathbf{a}_{M_2}$ be $K_1\times K_2$ filters supported on ${\mathbb K}$, and we introduce
\begin{align}
{\boldsymbol \mW}&=\left[{\boldsymbol \mS}_{\mathbf{a}_1[-\cdot]}^T,{\boldsymbol \mS}_{\mathbf{a}_2[-\cdot]}^T,\cdots,{\boldsymbol \mS}_{\mathbf{a}_{M_2}[-\cdot]}^T\right]^T,\label{OurAnalysis}\\
{\boldsymbol \mW}^*&=\left[{\boldsymbol \mS}_{\overline{\mathbf{a}}_1},{\boldsymbol \mS}_{\overline{\mathbf{a}}_2},\cdots,{\boldsymbol \mS}_{\overline{\mathbf{a}}_{M_2}}\right],\label{OurSynthesis}
\end{align}
where ${\boldsymbol \mS}_{\mathbf{a}}$ is a discrete convolution under the periodic boundary condition defined as \cref{DiscreteConv}. In other words, both ${\boldsymbol \mW}$ and ${\boldsymbol \mW}^*$ are concatenations of discrete convolutions.

\begin{theorem}\label{Th1} Let $u(\mathbf{x})$ be defined as in \cref{uModel} and $\boldsymbol{v}={\mathscr F}(u)\big|_{L^{-1}{\mathbb O}}$. Let ${\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)\in{\Bbb C}^{2M_1\times M_2}$ be a two-fold Hankel matrix with ${\boldsymbol\Lambda}$ defined as \cref{Lambda}. Assume that $\mathrm{rank}\left({\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)\right)=r\ll2M_1\wedge M_2$. Considering its full SVD ${\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)={\boldsymbol X}\boldsymbol{\Sigma}{\boldsymbol Y}^*$, we let $\mathbf{a}_j=M_2^{-1/2}{\boldsymbol Y}^{(:,j)}$ by reformulating ${\boldsymbol Y}^{(:,j)}\in{\Bbb C}^{M_2}$ into a $K_1\times K_2$ filter supported on ${\mathbb K}$. Then ${\boldsymbol \mW}$ in \cref{OurAnalysis} defined by using these $\mathbf{a}_1,\cdots,\mathbf{a}_{M_2}$ satisfies
\begin{align}\label{LowRankHankelTightFrame}
{\boldsymbol \mW}^*{\boldsymbol \mW}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)=\sum_{j=1}^{M_2}{\boldsymbol \mS}_{\overline{\mathbf{a}}_j}\left({\boldsymbol \mS}_{\mathbf{a}_j[-\cdot]}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)\right)={\boldsymbol\Lambda}\boldsymbol{v},
\end{align}
and for $j=r+1,\cdots,M_2,$ we have
\begin{align}\label{TightFrameSparse}
\left({\boldsymbol \mS}_{\mathbf{a}_j[-\cdot]}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)\right)[{\boldsymbol k}]={\mathbf 0},~~~~~{\boldsymbol k}\in{\mathbb O}:{\mathbb K},
\end{align}
where the discrete convolution ${\boldsymbol \mS}_{\mathbf{a}}$ acts on ${\boldsymbol\Lambda}_1\boldsymbol{v}$ and ${\boldsymbol\Lambda}_2\boldsymbol{v}$ separately\footnote{In fact, \cref{TightFrameSparse} is independent from the boundary condition of the discrete convolution, as ${\mathbb O}:{\mathbb K}$ is the grid where the boundary condition is inactive.}. Consequently, if ${\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)$ is of low rank, then its right singular vectors construct a tight frame under which ${\boldsymbol\Lambda}\boldsymbol{v}$ is sparsely represented.
\end{theorem}

In words, \cref{Th1} implies that a tight frame under which the weighted k-space data ${\boldsymbol\Lambda}\boldsymbol{v}$ is sparse can be explicitly constructed from the SVD of the corresponding Hankel matrix, and the sparsity of the tight frame coefficients is grouped according to filters. Even though the Hankel matrix is not fully known in the CS-MRI setting since the k-space data is undersampled, \cref{Th1} does motivate us to learn a better tight frame from the undersampled and noisy k-space data such that ${\boldsymbol\Lambda}\boldsymbol{v}$ can be sparsely represented better. More precisely, we remove the group sparsity pattern in the canonical coefficients ${\boldsymbol \mW}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)$, and find the coefficients which are as sparse as possible. This leads to solving
\begin{align}\label{AnalysisVersion}
\begin{split}
&~~~\min_{\boldsymbol{v},{\boldsymbol \mW}}~\frac{1}{2}\left\|{\boldsymbol \mR}_{{\mathbb M}}\boldsymbol{v}-{\boldsymbol f}\right\|_2^2+\gamma\left\|{\boldsymbol \mW}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)\right\|_0\\
&\text{subject to}~~\boldsymbol{v}\in{\mathcal C}~~\text{and}~~{\boldsymbol \mW}^*{\boldsymbol \mW}={\boldsymbol \mI}.
\end{split}
\end{align}
Here, ${\boldsymbol \mW}$ is a tight frame transform defined as in \cref{OurAnalysis}, i.e. the concatenation of the convolutions with filter banks $\left\{\mathbf{a}_1,\mathbf{a}_2,\cdots,\mathbf{a}_{M_2}\right\}$ which have to be learned, and ${\mathcal C}$ is a constraint set imposing the boundedness of $\boldsymbol{v}$. By relaxing the sparsity of tight frame coefficients over the entire range of ${\boldsymbol \mW}$ (not necessarily in groups as in \cref{TightFrameSparse}), we expect to achieve more flexibility, thereby leading to the improvements in restoration performance.

However, the direct penalization on the canonical coefficients ${\boldsymbol \mW}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)$ requires an iterative solver, which will be too expensive in constructing an adaptive tight frame \cite{J.F.Cai2014}. Therefore, by introducing an auxiliary variable ${\boldsymbol c}=\left({\boldsymbol c}_1,{\boldsymbol c}_2\right)$, we approximate \cref{AnalysisVersion} by
\begin{align}\label{ProposedCSMRIModel}
\begin{split}
&~~~\min_{\boldsymbol{v},{\boldsymbol c},{\boldsymbol \mW}}~\frac{1}{2}\left\|{\boldsymbol \mR}_{{\mathbb M}}\boldsymbol{v}-{\boldsymbol f}\right\|_2^2+\frac{\mu}{2}\left\|{\boldsymbol \mW}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)-{\boldsymbol c}\right\|_2^2+\gamma\left\|{\boldsymbol c}\right\|_0\\
&\text{subject to}~~\boldsymbol{v}\in{\mathcal C}~~\text{and}~~{\boldsymbol \mW}^*{\boldsymbol \mW}={\boldsymbol \mI},
\end{split}
\end{align}
to learn the tight frame system and reconstruct the fully sampled k-space data $\boldsymbol{v}$ on the grid ${\mathbb O}$ simultaneously.

\begin{remark}\label{RK1} For the constraint set ${\mathcal C}$, we choose
\begin{align}\label{ConstraintSet}
{\mathcal C}=\left\{\boldsymbol{v}\in{\mathcal V}_2:|\boldsymbol{v}[{\boldsymbol k}]|\leq R~~~\text{for all}~~{\boldsymbol k}\in{\mathbb O}\right\}
\end{align}
with a sufficiently large $R>0$, and such a choice comes from the following: if $u\in L_1({\Bbb R}^2)$ is modeled as in \cref{uModel}, then the (distributional) gradient satisfies
\begin{align*}
\nabla u=\left(\partial_1 u,\partial_2 u\right)=-\sum_{j=1}^J\alpha_j\boldsymbol{n}_j\mathrm{d}\boldsymbol{\sigma}\big|_{\partial\Omega_j},
\end{align*}
where $\boldsymbol{n}_j$ denotes the outward normal vector on $\partial\Omega_j$, and $\boldsymbol{\sigma}$ denotes the surface measure. Since $\nabla u$ is a Radon vector measure on ${\Bbb R}^2$ supported on $\Gamma$, its Fourier transform is written as the Fourier transform of a measure (e.g. \cite{Folland1999}):
\begin{align*}
{\mathscr F}(\nabla u)({\boldsymbol{\xi}})=-\sum_{j=1}^J\alpha_j\int_{\partial\Omega_j}e^{-2\pi i{\boldsymbol{\xi}}\cdot\mathbf{x}}\boldsymbol{n}_j(\mathbf{x})\mathrm{d}\boldsymbol{\sigma}(\mathbf{x}),~~~~~{\boldsymbol{\xi}}\in{\Bbb R}^2,
\end{align*}
which leads to
\begin{align*}
\left|{\mathscr F}(\nabla u)({\boldsymbol{\xi}})\right|\leq\sum_{j=1}^J|\alpha_j|\boldsymbol{\sigma}(\Omega_j)<\infty~~~~~\text{for all}~~{\boldsymbol{\xi}}\in{\Bbb R}^2.
\end{align*}
Hence, $\widehat{u}({\boldsymbol{\xi}})$ decays at infinity in the sense that
\begin{align*}
\sup_{{\boldsymbol{\xi}}\in{\Bbb R}^2}\left|{\boldsymbol{\xi}}\right|\left|\widehat{u}({\boldsymbol{\xi}})\right|<\infty,
\end{align*}
and this implies that $\widehat{u}({\boldsymbol{\xi}})$ (and whence $\boldsymbol{v}[{\boldsymbol k}]$ since we want $\boldsymbol{v}={\mathscr F}(u)\big|_{L^{-1}{\mathbb O}}$) has to be bounded. In addition, as mentioned in e.g. \cite{C.Bao2016}, the boundedness constraint is required for the convergence guarantee of the proximal alternating scheme to solve \cref{ProposedCSMRIModel}. That being said, numerically, the constraint set ${\mathcal C}$ 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 ${\mathcal C}$.
\end{remark}

Since our data driven tight frame approach is inspired by the SVD of a low rank ${\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)={\boldsymbol X}\boldsymbol{\Sigma}{\boldsymbol Y}^*$, we expect the learned filters $\left\{\mathbf{a}_1,\mathbf{a}_2,\cdots,\mathbf{a}_{M_2}\right\}$ are related to the annihilating filters. More precisely, we expect ${\boldsymbol Y}=M_2^{1/2}\left[\begin{array}{cccc}
\mathbf{a}_1&\mathbf{a}_2&\cdots&\mathbf{a}_{M_2}
\end{array}\right]$, so that $M_2^{1/2}\mathbf{a}_{r+1},\cdots,M_2^{1/2}\mathbf{a}_{M_2}$ form an orthonormal basis for the null space of ${\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)$. Hence, following the arguments in \cite{G.Ongie2016}, we can expect that the edge sets can also be restored from the learned filter banks as well. To see this, we let $\varphi_j$ be the trigonometric polynomial defined as
\begin{align*}
\varphi_j(\mathbf{x})=L^{-1}\sum_{{\boldsymbol k}\in{\mathbb K}}\mathbf{a}_j[-{\boldsymbol k}]e^{2\pi i{\boldsymbol k}\cdot\mathbf{x}/L},~~~~~j=r+1,\cdots,M_2,
\end{align*}
where the filters are flipped because we use the two-fold Hankel matrix formulation. Given that the conditions in \cite[Theorem 3.4]{G.Ongie2016} are satisfied, we can write $\varphi_j=\eta_j\varphi$ where $\varphi$ is the minimal polynomial in \cref{MajorAssumption} and $\eta_j$ is another $L$-periodic trigonometric polynomial. Following the similar arguments to \cite[Proposition 5.3]{G.Ongie2016}, for $\mathbf{x}\in\left[-L/2,L/2\right)^2$, we have
\begin{align}\label{EdgeEstimate}
\left\|{\boldsymbol V}\bV^*{\bf e}_{\mathbf{x}}\right\|_2=\left(\sum_{j=r+1}^{M_2}\left|\varphi_j(\mathbf{x})\right|^2\right)^{1/2}=\left(\sum_{j=r+1}^{M_2}\left|\eta_j(\mathbf{x})\varphi(\mathbf{x})\right|^2\right)^{1/2},
\end{align}
where
\begin{align*}
{\boldsymbol V}=M_2^{1/2}\left[\begin{array}{ccc}
\mathbf{a}_{r+1}[-\cdot]&\cdots&\mathbf{a}_{M_2}[-\cdot]
\end{array}\right]~~~\text{and}~~~{\bf e}_{\mathbf{x}}[{\boldsymbol k}]=L^{-1}e^{2\pi i{\boldsymbol k}\cdot\mathbf{x}/L},~~~{\boldsymbol k}\in{\mathbb K}.
\end{align*}
Hence, \cref{EdgeEstimate} shows that the learned filter banks are indeed related to the annihilating polynomials in the image domain.

Finally, we would like to mention that the application of dictionary learning to the CS-MRI reconstruction is already explored in \cite{S.Ravishankar2015}. Specifically, the authors in \cite{S.Ravishankar2015} proposed the following (unitary) sparsifying transform learning based CS-MRI (TLMRI) reconstruction model
\begin{align}\label{TLMRI}
\begin{split}
&~~~\min_{\boldsymbol{u},{\boldsymbol W},{\boldsymbol C}}~\mu\left\|{\boldsymbol \mR}_{{\mathbb M}}{\boldsymbol \msF}\boldsymbol{u}-{\boldsymbol f}\right\|_2^2+\sum_{j=1}^{N_p}\left\|{\boldsymbol W}{\boldsymbol \mP}_j\boldsymbol{u}-{\boldsymbol C}^{(:,j)}\right\|_2^2\\
&\text{subject to}~~\left\|\boldsymbol{u}\right\|_2\leq R,~~\left\|{\boldsymbol C}\right\|_0\leq r~~\text{and}~~{\boldsymbol W}^*{\boldsymbol W}={\boldsymbol I},
\end{split}
\end{align}
where ${\boldsymbol \mP}_j$ represents an operator that extracts an $K\times K$ patch of $\boldsymbol{u}$ as a vector, $N_p$ denotes the number of patches, and ${\boldsymbol W}$ denotes a square sparsifying transform for the patches of $\boldsymbol{u}$. Indeed, it seems that our proposed model \cref{ProposedCSMRIModel} can also be classified into the dictionary learning based approach due to the data driven tight frame regularization. However, unlike the above TLMRI model \cref{TLMRI} which considers the sparse structure of (the patches of) the underlying discrete image, the sparse regularization via data driven tight frame in \cref{ProposedCSMRIModel} is a relaxation of the low rank two-fold Hankel matrix completion in the frequency domain, reflecting the piecewise constant function in a \emph{continuous domain} via annihilating filters. Hence, \cref{TLMRI} is an \emph{on-the-grid} approach, whereas \cref{ProposedCSMRIModel} is an \emph{off-the-grid} approach.

\subsection{Alternating minimization algorithm}\label{PAMAlgorithm}

To solve \cref{ProposedCSMRIModel}, we use the proximal alternating minimization (PAM) algorithm introduced in \cite{H.Attouch2010}. We initialize $\boldsymbol{v}_0={\boldsymbol \mP}_{{\mathcal C}}\left({\boldsymbol \mR}_{{\mathbb M}}{\boldsymbol f}\right)$ where
\begin{align}\label{ProjectOntoC}
{\boldsymbol \mP}_{{\mathcal C}}\left(\boldsymbol{v}\right)[{\boldsymbol k}]=\min\left\{\left|\boldsymbol{v}[{\boldsymbol k}]\right|,R\right\}\exp\left\{i\arg\boldsymbol{v}[{\boldsymbol k}]\right\}.
\end{align}
For ${\boldsymbol \mW}_0$, we can use the SVD of ${\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}_0\right)$. However, to reduce the computational burden, we use the following way. Noting that we in general densely undersample the k-space data on the low frequencies based on the calibrationless sampling strategy, we restrict ${\boldsymbol\Lambda}\boldsymbol{v}_0$ onto the smaller grid, for instance, ${\mathbb O}'=\left\{-\lfloor N/4\rfloor,\cdots,\lfloor N/4-1\rfloor\right\}^2$. Then we compute the SVD ${\boldsymbol \mH}\left({\boldsymbol \mR}_{{\mathbb O}'}{\boldsymbol\Lambda}\boldsymbol{v}_0\right)={\boldsymbol X}_0\boldsymbol{\Sigma}_0{\boldsymbol Y}_0^*$, and set ${\boldsymbol \mW}_0$ via
\begin{align}\label{InitialFilterBanks}
\mathbf{a}_{0,j}=M_2^{-1/2}{\boldsymbol Y}_0^{(:,j)}~~~~~j=1,\cdots,M_2.
\end{align}
Using the initial filter banks as in \cref{InitialFilterBanks}, we initialize ${\boldsymbol c}_0$ via
\begin{align}\label{Initialize}
{\boldsymbol c}_0=\left[\begin{array}{cccccc}
\left[{\boldsymbol \mS}_{\mathbf{a}_{0,1}[-\cdot]}\left({\boldsymbol\Lambda}\boldsymbol{v}_0\right)\right]^T&\cdots&\left[{\boldsymbol \mS}_{\mathbf{a}_{0,r}[-\cdot]}\left({\boldsymbol\Lambda}\boldsymbol{v}_0\right)\right]^T&{\mathbf 0}&\cdots&{\mathbf 0}
\end{array}\right]^T,
\end{align}
where $r$ is an estimated rank.  After the initializations, we optimize $(\boldsymbol{v},{\boldsymbol c},{\boldsymbol \mW})$ alternatively, as summarized in \cref{Alg1}.

\begin{algorithm}[tp!]
\caption{Proximal Alternating Minimization Algorithm for \cref{ProposedCSMRIModel}}\label{Alg1}
\begin{algorithmic}
\STATE{\textbf{Initialization:} $\boldsymbol{v}_0$, ${\boldsymbol c}_0$, ${\boldsymbol \mW}_0$}
\FOR{$n=0$, $1$, $2$, $\cdots$}
\STATE{\textbf{(1)} Optimize $\boldsymbol{v}$:
\begin{align}\label{vsubprob}
\boldsymbol{v}_{n+1}=\operatornamewithlimits{argmin}_{\boldsymbol{v}\in{\mathcal C}}~\frac{1}{2}\left\|{\boldsymbol \mR}_{{\mathbb M}}\boldsymbol{v}-{\boldsymbol f}\right\|_2^2+\frac{\mu}{2}\left\|{\boldsymbol \mW}_n\left({\boldsymbol\Lambda}\boldsymbol{v}\right)-{\boldsymbol c}_n\right\|_2^2+\frac{\beta_1}{2}\left\|\boldsymbol{v}-\boldsymbol{v}_n\right\|_2^2.
\end{align}
\textbf{(2)} Optimize ${\boldsymbol c}$:
\begin{align}\label{csubprob}
{\boldsymbol c}_{n+1}=\operatornamewithlimits{argmin}_{{\boldsymbol c}}~\gamma\left\|{\boldsymbol c}\right\|_0+\frac{\mu}{2}\left\|{\boldsymbol c}-{\boldsymbol \mW}_{n}\left({\boldsymbol\Lambda}\boldsymbol{v}_{n+1}\right)\right\|_2^2+\frac{\beta_2}{2}\left\|{\boldsymbol c}-{\boldsymbol c}_n\right\|_2^2.
\end{align}
\textbf{(3)} Optimize ${\boldsymbol \mW}$:
\begin{align}\label{Wsubprob}
{\boldsymbol \mW}_{n+1}=\operatornamewithlimits{argmin}_{{\boldsymbol \mW}^*{\boldsymbol \mW}={\boldsymbol \mI}}~\frac{\mu}{2}\left\|{\boldsymbol \mW}\left({\boldsymbol\Lambda}\boldsymbol{v}_{n+1}\right)-{\boldsymbol c}_{n+1}\right\|_2^2+\frac{\beta_3}{2}\left\|{\boldsymbol \mW}-{\boldsymbol \mW}_n\right\|_F^2.
\end{align}}
\ENDFOR
\end{algorithmic}
\end{algorithm}

It is easy to see that each subproblem in \cref{Alg1} has a closed form solution. Since ${\boldsymbol \mW}_n^*{\boldsymbol \mW}_n={\boldsymbol \mI}$ for all $n$, the solution to \cref{vsubprob} is given by
\begin{align}\label{vsubexplicit}
\begin{split}
\boldsymbol{v}_{n+1/2}&=\left({\boldsymbol \mR}_{{\mathbb M}}^T{\boldsymbol \mR}_{{\mathbb M}}+\mu{\boldsymbol\Lambda}^*{\boldsymbol\Lambda}+\beta_1{\boldsymbol \mI}\right)^{-1}\left({\boldsymbol \mR}_{{\mathbb M}}^T{\boldsymbol f}+\mu{\boldsymbol\Lambda}^*{\boldsymbol \mW}_n^*{\boldsymbol c}_n+\beta_1\boldsymbol{v}_n\right),\\
\boldsymbol{v}_{n+1}&={\boldsymbol \mP}_{{\mathcal C}}\left(\boldsymbol{v}_{n+1/2}\right),
\end{split}
\end{align}
where ${\boldsymbol \mR}_{{\mathbb M}}^T$ is the adjoint of ${\boldsymbol \mR}_{{\mathbb M}}$ (i.e. the zero-filling operator), and ${\boldsymbol \mP}_{{\mathcal C}}$ is defined as \cref{ProjectOntoC}. It is worth noting that since ${\boldsymbol \mR}_{{\mathbb M}}+\mu{\boldsymbol\Lambda}^*{\boldsymbol\Lambda}+\beta_1{\boldsymbol \mI}$ is a diagonal matrix, no matrix inversion is needed.

For the remaining subproblems \cref{csubprob,Wsubprob}, we note that while the coefficients in \cref{DDTFModel} are intermediate variables to construct tight frame filter banks, the coefficients in \cref{ProposedCSMRIModel} are used to update $\boldsymbol{v}$. Hence, we solve \cref{csubprob,Wsubprob} in the following way. The solution to \cref{csubprob} is expressed as
\begin{align}\label{csubexplicit}
{\boldsymbol c}_{n+1}={\boldsymbol \mT}_{\sqrt{2\gamma/(\mu+\beta_2)}}\left(\frac{\mu{\boldsymbol \mW}_n\left({\boldsymbol\Lambda}\boldsymbol{v}_{n+1}\right)+\beta_2{\boldsymbol c}_n}{\mu+\beta_2}\right),
\end{align}
where the hard thresholding ${\boldsymbol \mT}$ for ${\boldsymbol c}=\left({\boldsymbol c}_1,{\boldsymbol c}_2\right)$ is defined as
\begin{align}\label{HardThresholding}
{\boldsymbol \mT}_{\gamma}\left({\boldsymbol c}\right)=\left({\boldsymbol \mT}_{\gamma}\left({\boldsymbol c}_1\right),{\boldsymbol \mT}_{\gamma}\left({\boldsymbol c}_2\right)\right),~~~\text{where}~~~{\boldsymbol \mT}_{\gamma}\left({\boldsymbol c}_j\right)[{\boldsymbol k}]=\left\{\begin{array}{cl}
{\boldsymbol c}_j[{\boldsymbol k}]~&\text{if}~|{\boldsymbol c}_j[{\boldsymbol k}]|>\lambda\\
0~&\text{otherwise.}
\end{array}\right.
\end{align}
To solve \cref{Wsubprob}, let ${\boldsymbol H}_{n+1}={\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}_{n+1}\right)\in{\Bbb C}^{2M_1\times M_2}$ and let ${\boldsymbol A}\in{\Bbb C}^{M_2\times M_2}$ be a matrix whose column vectors are concatenations of the filters $\mathbf{a}_1,\cdots,\mathbf{a}_{M_2}$. Then we have
\begin{align}
{\boldsymbol A}_{n+1}&=\operatornamewithlimits{argmin}_{{\boldsymbol A}\bsA^*=M_2^{-1}{\boldsymbol I}}~\frac{\mu}{2}\left\|{\boldsymbol H}_{n+1}{\boldsymbol A}-{\boldsymbol C}_{n+1}\right\|_F^2+\frac{\beta_3}{2}\left\|{\boldsymbol A}-{\boldsymbol A}_n\right\|_F^2,\label{Dupdate}
\end{align}
where ${\boldsymbol C}=\left[\begin{array}{cc}
{\boldsymbol C}_1^T&{\boldsymbol C}_2^T
\end{array}\right]^T\in{\Bbb C}^{2M_1\times M_2}$, and each ${\boldsymbol C}_j\in{\Bbb C}^{M_1\times M_2}$ is the corresponding matrix formulation of frame coefficients ${\boldsymbol c}_j$. Then the closed form solution to \cref{Dupdate} is expressed as
\begin{align}\label{Wsubexplicit}
{\boldsymbol A}_{n+1}=M_2^{-1/2}{\boldsymbol X}_n{\boldsymbol Y}_n^*~~\text{with}~~{\boldsymbol X}_n\boldsymbol{\Sigma}_n{\boldsymbol Y}_n^*={\boldsymbol H}_{n+1}^*{\boldsymbol C}_{n+1}+\mu^{-1}\beta_3{\boldsymbol A}_n~~\text{is the SVD}.
\end{align}
Therefore, we only compute the SVD of ${\boldsymbol H}_{n+1}^*{\boldsymbol C}_{n+1}+\mu^{-1}\beta_3{\boldsymbol A}_n\in{\Bbb C}^{M_2\times M_2}$, leading to the computational efficiency over directly minimizing the rank of ${\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)\in{\Bbb C}^{2M_1\times M_2}$.

Note that, based on the frameworks in \cite{H.Attouch2010,C.Bao2016,J.Bolte2014,Kurdyk1998,Lojasiewicz1995,Y.Xu2013}, it can be verified that \cref{Alg1} globally converges to the critical point, as the minimization of a real-valued complex variable objective function in \cref{ProposedCSMRIModel} is based on the Wirtinger calculus (i.e. by identifying ${\Bbb C}\simeq{\Bbb R}^2$).

For the computational cost, we first mention that we do not explicitly formulate ${\boldsymbol H}_{n+1}={\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}_{n+1}\right)$ at each iteration step. Without loss of generality, we assume the square patch, i.e. $K_1=K_2=K$, $M_1=(N-K+1)^2$ and $M_2=K^2$. Then ${\boldsymbol \mW}_n\left({\boldsymbol\Lambda}\boldsymbol{v}_{n+1}\right)$ and ${\boldsymbol H}_{n+1}^*{\boldsymbol C}_{n+1}$ are computed by using $4K^2$ ($2K^2$ for ${\boldsymbol \mW}_n\left({\boldsymbol\Lambda}\boldsymbol{v}_{n+1}\right)$ and $2K^2$ for ${\boldsymbol H}_{n+1}^*{\boldsymbol C}_{n+1}$) two dimensional convolutions/fast Hankel matrix-vector multiplications \cite{L.Lu2015} directly from ${\boldsymbol\Lambda}\boldsymbol{v}_{n+1}$, requiring $O(K^2N^2\log N)$ operations. In addition, since the hard thresholding ${\boldsymbol \mT}$ is an elementwise operation, we can update ${\boldsymbol H}_{n+1}^*{\boldsymbol C}_{n+1}^{(:,j)}+\mu^{-1}\beta_3{\boldsymbol A}_n^{(:,j)}$ directly after updating the $j$th element of ${\boldsymbol c}_{n+1}$. Hence, together with the SVD of ${\boldsymbol H}_{n+1}^*{\boldsymbol C}_{n+1}+\mu^{-1}\beta_3{\boldsymbol A}_n\in{\Bbb C}^{K^2\times K^2}$ requiring $O(K^6)$ operations, we need $O(K^2N^2\log N+K^2N^2+K^6)$ operations for \cref{csubexplicit,Wsubexplicit}. For \cref{vsubexplicit}, note that ${\boldsymbol \mW}_n^*{\boldsymbol c}_n$ can be computed by using $2K^2$ two dimensional convolutions, requiring $O(K^2N^2\log N)$ operations. Then since the remaining operations are all pointwise multiplications and additions, \cref{vsubexplicit} requires $O(K^2N^2\log N+N^2)$ operations. Combining these all together, the computational cost of \cref{Alg1} is $O(K^2N^2\log N+K^2N^2+K^6)$.

We further mention that our PAM algorithm seems to be similar to the block coordinate descent (BCD) algorithm presented in \cite{S.Ravishankar2015} to solve the TLMRI model \cref{TLMRI}. In fact, when we set $\beta_1=\beta_2=\beta_3=0$, \cref{Alg1} indeed reduces to the BCD algorithm. In addition, ${\boldsymbol W}$ subproblem in \cref{TLMRI} can be solved similarly to \cref{Wsubexplicit}. Nevertheless, there are differences between \cref{Alg1} and the BCD algorithm in \cite{S.Ravishankar2015}. First of all, since \cref{ProposedCSMRIModel} aims to restore fully sampled k-space data, we do not explicitly use fast Fourier transform (FFT) during the iteration, unlike \cref{TLMRI} which requires to use FFT several times in updating the image $\boldsymbol{u}$. Second, as previously mentioned, it is not necessary to explicitly formulate the two-fold Hankel matrix for \cref{csubexplicit,Wsubexplicit}, while ${\boldsymbol C}$ subproblem and ${\boldsymbol W}$ subproblem in \cref{TLMRI}, i.e. the so-called \emph{sparse coding} subproblem requires to extract each patch of the image and store it as a column vector. Finally, even though the subsequence convergence of the BCD algorithm for \cref{TLMRI} is established in the works including \cite{S.Ravishankar2015,B.Wen2017}, this does not necessarily imply the global convergence of the iterates. Meanwhile, as previously mentioned, \cref{Alg1} for \cref{ProposedCSMRIModel} guarantees the global convergence of the iterates. Even though the global convergence is not crucial to the applications in image restoration in general, we propose the PAM algorithm with a global convergence guarantee, to balance both theoretical interest and the needs from applications.

\section{Experimental results}\label{Experiments}

In this section, we present the experimental results on the phantom image and the real MR image used in \cite{G.Ongie2016}, to compare the proposed DDTF based CS-MRI model \cref{ProposedCSMRIModel} with several existing methods. We choose to compare with the total variation (TV) model \cite{M.Lustig2007}
\begin{align}\label{TVModel}
\min_{\boldsymbol{u}}~\frac{1}{2}\left\|{\boldsymbol \mR}_{{\mathbb M}}{\boldsymbol \msF}\boldsymbol{u}-{\boldsymbol f}\right\|_2^2+\gamma\left\|\boldsymbol \nabla\boldsymbol{u}\right\|_1,
\end{align}
the Haar framelet (Haar) model (e.g. \cite{J.F.Cai2012})
\begin{align}\label{HaarModel}
\min_{\boldsymbol{u}}~\frac{1}{2}\left\|{\boldsymbol \mR}_{{\mathbb M}}{\boldsymbol \msF}\boldsymbol{u}-{\boldsymbol f}\right\|_2^2+\gamma\left\|{\boldsymbol \mW}\boldsymbol{u}\right\|_1,
\end{align}
and the TLMRI model \cref{TLMRI} where a learned transform ${\boldsymbol W}$ is a unitary flipping and rotation invariant sparsifying transform \cite{B.Wen2017}. Both \cref{TVModel,HaarModel} are solved by the split Bregman method \cite{J.F.Cai2009/10,T.Goldstein2009}, and \cref{TLMRI} is solved by the BCD algorithm in \cite{S.Ravishankar2015,B.Wen2017}. We also compare with the following Schatten $p$-norm minimization model
\begin{align}\label{Schatten}
\min_{\boldsymbol{v}\in{\mathcal C}}~\left\|{\boldsymbol \mR}_{{\mathbb M}}\boldsymbol{v}-{\boldsymbol f}\right\|_2^2+\gamma\left\|{\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)\right\|_p^p,
\end{align}
with $p=0$, $0.5$, and $1$, solved by the generic iterative reweighted annihilating filters (GIRAF) method \cite{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.

\begin{table}[tp!]
\centering
\caption{Parameter selection for each dataset. For TV \cref{TVModel}, Haar \cref{HaarModel}, and GIRAF \cref{Schatten}, $\mu$ refers to the internal parameter for the split Bregman algorithm.}\label{TableParameter}
\vspace{-0.20cm}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
\hline
Dataset&Model&$K$&$r$&$\mu$&$\gamma$&$\beta_1$&$\beta_2$&$\beta_3$\\ \hline\hline
\multirow{5}{*}{Phantom}&TV \cref{TVModel}&$\cdot$&$\cdot$&$10$&$0.05$&$\cdot$&$\cdot$&$\cdot$\\ \cline{2-9}
&Haar \cref{HaarModel}&$\cdot$&$\cdot$&$10$&$0.1$&$\cdot$&$\cdot$&$\cdot$\\ \cline{2-9}
&TLMRI \cref{TLMRI}&$8$&$198403$&$15.2588$&$\cdot$&$\cdot$&$\cdot$&$\cdot$\\ \cline{2-9}
&GIRAF \cref{Schatten}&$25$&$\cdot$&$4$&$10^{-6}$&$\cdot$&$\cdot$&$\cdot$\\ \cline{2-9}
&DDTF \cref{ProposedCSMRIModel}&$25$&$500$&$0.1$&$10$&$10^{-4}$&$10^{-4}$&$10^{-4}$\\ \hline\hline
\multirow{5}*{Real MR}&TV \cref{TVModel}&$\cdot$&$\cdot$&$10$&$0.01$&$\cdot$&$\cdot$&$\cdot$\\ \cline{2-9}
&Haar \cref{HaarModel}&$\cdot$&$\cdot$&$10$&$0.025$&$\cdot$&$\cdot$&$\cdot$\\ \cline{2-9}
&TLMRI \cref{TLMRI}&$8$&$396806$&$15.2588$&$\cdot$&$\cdot$&$\cdot$&$\cdot$\\ \cline{2-9}
&GIRAF \cref{Schatten}&$45$&$\cdot$&$2$&$10^{-6}$&$\cdot$&$\cdot$&$\cdot$\\ \cline{2-9}
&DDTF \cref{ProposedCSMRIModel}&$45$&$1620$&$0.05$&$5$&$10^{-4}$&$10^{-4}$&$10^{-4}$\\ \hline
\end{tabular}
\end{table}

In \cref{ProposedCSMRIModel,Schatten}, we choose $R=|{\boldsymbol f}[{\mathbf 0}]|$ if ${\mathbf 0}\in{\mathbb M}$ and $R=10^8$ otherwise for the constraint set ${\mathcal C}$ in \cref{ConstraintSet}. In \cref{TLMRI}, we choose $R=10^5$. We use the forward difference for the discrete gradient $\boldsymbol \nabla$ in \cref{TVModel}, and ${\boldsymbol \mW}$ in \cref{HaarModel} is chosen to be the undecimated tensor product Haar framelet transform with $1$ level of decomposition \cite{B.Dong2013}. Both \cref{ProposedCSMRIModel,Schatten} use the $K\times K$ square patch to generate the two-fold Hankel matrix for simplicity. We choose $K$ to be no more than $0.2N$, based on \cite{G.Ongie2018}, and we choose $r\approx0.8K^2$ for the initialization \cref{Initialize} of \cref{ProposedCSMRIModel}, both of which depend on the geometry of the target image. For \cref{TLMRI}, we use $8\times 8$ patches. The detailed choice of the remaining regularization parameters are summarized in \cref{TableParameter}. Empirically, we observe that $\mu\approx0.01\gamma$ is an appropriate choice for \cref{ProposedCSMRIModel}, and we observe that when $\mu$ is large, the restored k-space data $\boldsymbol{v}$ 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  \cref{TVModel,HaarModel,TLMRI}, the stopping criterion is
\begin{align}\label{StoppingOntheGrid}
\frac{\left\|\boldsymbol{u}_{n+1}-\boldsymbol{u}_n\right\|_2}{\left\|\boldsymbol{u}_n\right\|_2}\leq\varepsilon:=2\times10^{-4},
\end{align}
and for the off-the-grid approaches \cref{ProposedCSMRIModel,Schatten}, we use
\begin{align}\label{StoppingOfftheGrid}
\frac{\left\|\boldsymbol{v}_{n+1}-\boldsymbol{v}_n\right\|_2}{\left\|\boldsymbol{v}_n\right\|_2}\leq\varepsilon.
\end{align}
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), the high frequency error norm (HFEN) \cite{S.Ravishankar2011}, the number of iterations (NIters) to reach the stopping criterion, and the CPU time. Note that for the off-the-grid approaches \cref{ProposedCSMRIModel,Schatten}, the restored image is computed via the inverse DFT of the restored k-space data; see \cref{CSMRIComparisons}.

\subsection{Phantom experiments}\label{PhantomExperiments}

For the piecewise constant phantom experiments, we first compute the fully sampled k-space data using \cref{MRForward} with a piecewise constant function model \cref{uModel}. More precisely, given that $\Omega_j$'s are ellipses, \cref{uModel} can be explicitly written as
\begin{align}\label{uModelExp}
u(\mathbf{x})=\sum_{j=1}^J\alpha_j1_{B({\mathbf 0},1)}\left({\boldsymbol D}_j{\boldsymbol Q}_j\left(\mathbf{x}-\mathbf{y}_j\right)\right),
\end{align}
where ${\boldsymbol D}_j$'s are diagonal matrices, ${\boldsymbol Q}_j$'s are rotation matrices, and $B({\mathbf 0},1)$ is the unit disk. Following \cite{M.Guerquin-Kern2012}, we can compute the Fourier transform of \cref{uModelExp} by noting that
\begin{align*}
{\mathscr F}(1_{B({\mathbf 0},1)})({\boldsymbol{\xi}})=\frac{J_1(2\pi|{\boldsymbol{\xi}}|)}{|{\boldsymbol{\xi}}|}~~~\text{and}~~~{\mathscr F}(1_{B({\mathbf 0},1)})({\mathbf 0})=\pi,
\end{align*}
where $J_1$ is the first kind Bessel's function of order $1$. Then using the variable density random sampling method in \cite{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 \cref{PhantomDataSet}).

\begin{figure}[tp!]
\centering
\hspace{-0.1cm}\subfloat[Fully sampled]{\label{PhantomOriginalk}\includegraphics[width=3.0cm]{PhantomOriginalk.pdf}}\hspace{0.005cm}
\subfloat[Ground truth]{\label{PhantomOriginal}\includegraphics[width=3.0cm]{PhantomOriginal.pdf}}\hspace{0.005cm}
\subfloat[Sample mask]{\label{PhantomMask}\includegraphics[width=3.0cm]{PhantomMask.pdf}}\hspace{0.005cm}
\subfloat[Undersampled]{\label{PhantomUndersamplek2}\includegraphics[width=3.0cm]{PhantomUndersamplek.pdf}}\vspace{-0.20cm}
\caption{Dataset for the phantom experiments. Fully sampled k-space data, its inverse DFT as a ground truth, the undersampling mask, and the undersampled k-space data.}\label{PhantomDataSet}
\end{figure}

\begin{figure}[tp!]
\centering
\hspace{-0.1cm}\subfloat[$r=550$]{\label{PhantomDDTFEdge551}\includegraphics[width=3.5cm]{PhantomDDTFEdge551Rev.pdf}}\hspace{0.005cm}
\subfloat[$r=450$]{\label{PhantomDDTFEdge451}\includegraphics[width=3.5cm]{PhantomDDTFEdge451Rev.pdf}}\hspace{0.005cm}
\subfloat[$r=350$]{\label{PhantomDDTFEdge351}\includegraphics[width=3.5cm]{PhantomDDTFEdge351Rev.pdf}}\hspace{0.005cm}
\subfloat[$r=336$]{\label{PhantomDDTFEdge337}\includegraphics[width=3.5cm]{PhantomDDTFEdge337Rev.pdf}}\vspace{-0.20cm}
\caption{Images of edges estimated by using \cref{EdgeEstimate} with the filter banks in \cref{ProposedCSMRIModel} and various values of $r$. All estimated images are displayed in the window level $[0,1]$.}\label{PhantomEdges}
\end{figure}

\begin{table}[tp!]
\centering
\caption{Comparison of signal-to-noise ratio, high frequency error norm, iteration numbers, and CPU time for the phantom experiments. MaxIter denotes that the iteration reached the maximum number without satisfying the stopping criterion.}\label{PhantomTable}
\vspace{-0.20cm}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
\hline
\multirow{2}{*}{Indices}&\multirow{2}{*}{Zero fill}&\multirow{2}{*}{TV \cref{TVModel}}&\multirow{2}{*}{Haar \cref{HaarModel}}&\multirow{2}{*}{TLMRI \cref{TLMRI}}&\multicolumn{3}{|c|}{GIRAF \cref{Schatten}}&\multirow{2}{*}{DDTF \cref{ProposedCSMRIModel}}\\ \cline{6-8}
&&&&&$p=0$&$p=0.5$&$p=1$&\\ \hline\hline
SNR&$9.23$&$20.29$&$20.99$&$21.62$&$21.12$&$20.88$&$17.13$&$\textbf{26.66}$\\ \hline
HFEN&$0.5490$&$0.0992$&$0.0904$&$0.0675$&$0.0948$&$0.0999$&$0.1828$&$\textbf{0.0572}$\\ \hline
NIters&$\cdot$&$210$&$244$&MaxIter&$18$&$20$&$15$&$154$\\ \hline
Time($\mathrm{s}$)&$\cdot$&$2.4$&$13.7$&$32612.9$&$15.0$&$16.2$&$12.2$&$1438.2$\\ \hline
\end{tabular}
\end{table}

\cref{PhantomEdges} displays the edges estimated by the proposed model \cref{ProposedCSMRIModel}. More precisely, using the obtained filter banks $\mathbf{a}_1,\cdots,\mathbf{a}_{K^2}$ at convergence, we compute $\left\|{\boldsymbol V}\bV^*{\bf e}_{\mathbf{x}}\right\|_2$ defined as \cref{EdgeEstimate}, by varying the rank cutoff $r=550, 450, 350, 336$, respectively. Here, $r=336$ is the ground truth rank cutoff described in \cite{G.Ongie2018}. At first glance, we can see that as the rank cutoff $r$ in \cref{EdgeEstimate} decreases, the image of $\left\|{\boldsymbol V}\bV^*{\bf e}_{\mathbf{x}}\right\|_2$ contains less artifacts. Most importantly, we can further observe that in any case, $\left\|{\boldsymbol V}\bV^*{\bf e}_{\mathbf{x}}\right\|_2$ is close to $0$ near the edges. This demonstrates that the data driven tight frame filter banks used in the proposed model \cref{ProposedCSMRIModel} are indeed related to the annihilating filter, as they are also able to estimate the edges of the piecewise constant function via \cref{EdgeEstimate}. Therefore, we conclude that the proposed sparse regularization via data driven tight frame provides another relaxation of the low rank two-fold Hankel matrix completion based on the annihilating filter.

\cref{PhantomTable} summarizes the SNR and the HFEN of the aforementioned restoration models, and \cref{PhantomResults} displays the visual comparisons with the zoom-in views in \cref{PhantomResultsZoom} and the error maps in \cref{PhantomResultsErrorMap}, respectively. We can see that the proposed data driven tight frame model \cref{ProposedCSMRIModel} consistently outperforms both the on-the-grid approaches (\cref{TVModel,HaarModel,TLMRI}) and the existing off-the-grid approaches \cref{Schatten} with a smaller error map. Noting that \cref{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 \cref{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 (\cref{PhantomTVError,PhantomHaarError}) compared to the off-the-grid approaches. Even though \cref{TLMRI} is able to restore the three small ellipses, the restored image is still overly smoothed compared to the proposed model with errors concentrated near the edges, as we can see from \cref{PhantomFRISTError,PhantomFRISTZoom}.

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 $\mathrm{rank}\left({\boldsymbol \mH}\left({\boldsymbol\Lambda}\boldsymbol{v}\right)\right)$ even in the fully sampled case as the weight matrix ${\boldsymbol\Lambda}$ 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 \cref{PhantomGIRAF0,PhantomGIRAFHalf,PhantomGIRAF1,PhantomGIRAF0Zoom,PhantomGIRAFHalfZoom,PhantomGIRAF1Zoom,PhantomGIRAF0Error,PhantomGIRAFHalfError,PhantomGIRAF1Error}. In contrast, the sparse approximation of ${\boldsymbol\Lambda}\boldsymbol{v}$ 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 \cref{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 \cref{PhantomTVk,PhantomHaark,PhantomFRISTk} that the restored k-space data by \cref{TVModel,HaarModel,TLMRI} 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 \cref{PhantomGIRAF0k,PhantomGIRAFHalfk,PhantomGIRAF1k}. In contrast, the proposed model \cref{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.

Finally, we note that even though \cref{ProposedCSMRIModel} requires more iterations than \cref{Schatten} for convergence, the number of iterations is less than \cref{TVModel,HaarModel,TLMRI}. We also note that the BCD algorithm of \cref{TLMRI} fails to meet the stopping criterion. It is likely that due to the subsequence convergence, the BCD algorithm for \cref{TLMRI} behaves inferior to other models, whose algorithms guarantee the global convergence. Moreover, it takes \cref{TLMRI} approximately $9\mathrm{hrs}$ to perform $600$ iterations. It is most likely that the generation of patch and the projection of a coefficient matrix onto the $\ell_0$ ball are most time consuming, affecting the CPU time. Meanwhile, despite the superior restoration performance, it takes \cref{ProposedCSMRIModel} approximately $24\mathrm{mins}$ to meet the stopping criterion due to a number of discrete convolutions. Even though this is shorter than \cref{TLMRI}, the proposed approach is still time-consuming compared to \cref{TVModel,HaarModel,Schatten}. Nevertheless, our approach is able to gain on the quality of the restored image.

\begin{figure}[tp!]
\centering
\hspace{-0.1cm}\subfloat[Fully sampled]{\label{PhantomOriginal2}\includegraphics[width=3.00cm]{PhantomOriginal.pdf}}\hspace{0.005cm}
\subfloat[Zero fill]{\label{PhantomZeroPad}\includegraphics[width=3.00cm]{PhantomZeroPad.pdf}}\hspace{0.005cm}
\subfloat[TV \cref{TVModel}]{\label{PhantomTV}\includegraphics[width=3.00cm]{PhantomTV.pdf}}\hspace{0.005cm}
\subfloat[Haar \cref{HaarModel}]{\label{PhantomHaar}\includegraphics[width=3.00cm]{PhantomHaar.pdf}}\hspace{0.005cm}
\subfloat[TLMRI \cref{TLMRI}]{\label{PhantomFRIST}\includegraphics[width=3.00cm]{PhantomFRIST.pdf}}\vspace{-0.20cm}\\
\subfloat[GIRAF$0$]{\label{PhantomGIRAF0}\includegraphics[width=3.00cm]{PhantomGIRAF0.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$0.5$]{\label{PhantomGIRAFHalf}\includegraphics[width=3.00cm]{PhantomGIRAFHalf.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$1$]{\label{PhantomGIRAF1}\includegraphics[width=3.00cm]{PhantomGIRAF1.pdf}}\hspace{0.005cm}
\subfloat[DDTF \cref{ProposedCSMRIModel}]{\label{PhantomDDTF}\includegraphics[width=3.00cm]{PhantomDDTF.pdf}}\vspace{-0.20cm}
\caption{Visual comparisons of each restoration models for the phantom experiments. All restored images are displayed in the window level $[0,1]$ for fair comparisons.}\label{PhantomResults}
\end{figure}

\begin{figure}[tp!]
\centering
\hspace{-0.1cm}\subfloat[Sample mask]{\label{PhantomMask2}\includegraphics[width=3.00cm]{PhantomMask.pdf}}\hspace{0.005cm}
\subfloat[Zero fill]{\label{PhantomZeroPadError}\includegraphics[width=3.00cm]{PhantomZeroPadError.pdf}}\hspace{0.005cm}
\subfloat[TV \cref{TVModel}]{\label{PhantomTVError}\includegraphics[width=3.00cm]{PhantomTVError.pdf}}\hspace{0.005cm}
\subfloat[Haar \cref{HaarModel}]{\label{PhantomHaarError}\includegraphics[width=3.00cm]{PhantomHaarError.pdf}}\hspace{0.005cm}
\subfloat[TLMRI \cref{TLMRI}]{\label{PhantomFRISTError}\includegraphics[width=3.00cm]{PhantomFRISTError.pdf}}\vspace{-0.20cm}\\
\subfloat[GIRAF$0$]{\label{PhantomGIRAF0Error}\includegraphics[width=3.00cm]{PhantomGIRAF0Error.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$0.5$]{\label{PhantomGIRAFHalfError}\includegraphics[width=3.00cm]{PhantomGIRAFHalfError.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$1$]{\label{PhantomGIRAF1Error}\includegraphics[width=3.00cm]{PhantomGIRAF1Error.pdf}}\hspace{0.005cm}
\subfloat[DDTF \cref{ProposedCSMRIModel}]{\label{PhantomDDTFError}\includegraphics[width=3.00cm]{PhantomDDTFError.pdf}}\vspace{-0.20cm}
\caption{Comparisons of error maps for the phantom experiments. All error maps are displayed in the window level $[0,0.2]$ for fair comparisons.}\label{PhantomResultsErrorMap}
\end{figure}

\begin{figure}[tp!]
\centering
\hspace{-0.1cm}\subfloat[Fully sampled]{\label{PhantomOriginalk2}\includegraphics[width=3.00cm]{PhantomOriginalk.pdf}}\hspace{0.005cm}
\subfloat[Zero fill]{\label{PhantomUndersamplek}\includegraphics[width=3.00cm]{PhantomUndersamplek.pdf}}\hspace{0.005cm}
\subfloat[TV \cref{TVModel}]{\label{PhantomTVk}\includegraphics[width=3.00cm]{PhantomTVk.pdf}}\hspace{0.005cm}
\subfloat[Haar \cref{HaarModel}]{\label{PhantomHaark}\includegraphics[width=3.00cm]{PhantomHaark.pdf}}\hspace{0.005cm}
\subfloat[TLMRI \cref{TLMRI}]{\label{PhantomFRISTk}\includegraphics[width=3.00cm]{PhantomFRISTk.pdf}}\vspace{-0.20cm}\\
\subfloat[GIRAF$0$]{\label{PhantomGIRAF0k}\includegraphics[width=3.00cm]{PhantomGIRAF0k.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$0.5$]{\label{PhantomGIRAFHalfk}\includegraphics[width=3.00cm]{PhantomGIRAFHalfk.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$1$]{\label{PhantomGIRAF1k}\includegraphics[width=3.00cm]{PhantomGIRAF1k.pdf}}\hspace{0.005cm}
\subfloat[DDTF \cref{ProposedCSMRIModel}]{\label{PhantomDDTFk}\includegraphics[width=3.00cm]{PhantomDDTFk.pdf}}\vspace{-0.20cm}
\caption{Comparisons of k-space data for the phantom experiments in the log scale. All restored k-space data are displayed in the window level $[0,9]$ for fair comparisons.}\label{PhantomResultsk}
\end{figure}

\begin{figure}[tp!]
\centering
\hspace{-0.1cm}\subfloat[Fully sampled]{\label{PhantomOriginalZoom}\includegraphics[width=3.00cm]{PhantomOriginalZoom.pdf}}\hspace{0.005cm}
\subfloat[Zero fill]{\label{PhantomZeroPadZoom}\includegraphics[width=3.00cm]{PhantomZeroPadZoom.pdf}}\hspace{0.005cm}
\subfloat[TV \cref{TVModel}]{\label{PhantomTVZoom}\includegraphics[width=3.00cm]{PhantomTVZoom.pdf}}\hspace{0.005cm}
\subfloat[Haar \cref{HaarModel}]{\label{PhantomHaarZoom}\includegraphics[width=3.00cm]{PhantomHaarZoom.pdf}}\hspace{0.005cm}
\subfloat[TLMRI \cref{TLMRI}]{\label{PhantomFRISTZoom}\includegraphics[width=3.00cm]{PhantomFRISTZoom.pdf}}\vspace{-0.20cm}\\
\subfloat[GIRAF$0$]{\label{PhantomGIRAF0Zoom}\includegraphics[width=3.00cm]{PhantomGIRAF0Zoom.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$0.5$]{\label{PhantomGIRAFHalfZoom}\includegraphics[width=3.00cm]{PhantomGIRAFHalfZoom.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$1$]{\label{PhantomGIRAF1Zoom}\includegraphics[width=3.00cm]{PhantomGIRAF1Zoom.pdf}}\hspace{0.005cm}
\subfloat[DDTF \cref{ProposedCSMRIModel}]{\label{PhantomDDTFZoom}\includegraphics[width=3.00cm]{PhantomDDTFZoom.pdf}}\vspace{-0.20cm}
\caption{Zoom-in views of \cref{PhantomResults}. The yellow arrows indicate the regions worth noticing.}\label{PhantomResultsZoom}
\end{figure}

\subsection{Real MR image experiments}\label{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 \cite{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 \cite{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 \cref{RealDataSet}.

\begin{figure}[tp!]
\centering
\hspace{-0.1cm}\subfloat[Fully sampled]{\label{RealOriginalk}\includegraphics[width=3.0cm]{RealOriginalk.pdf}}\hspace{0.005cm}
\subfloat[Ground truth]{\label{RealOriginal}\includegraphics[width=3.0cm]{RealOriginal.pdf}}\hspace{0.005cm}
\subfloat[Sample mask]{\label{RealMask}\includegraphics[width=3.0cm]{RealMask.pdf}}\hspace{0.005cm}
\subfloat[Undersampled]{\label{RealUndersamplek2}\includegraphics[width=3.0cm]{RealUndersamplek.pdf}}\vspace{-0.20cm}
\caption{Dataset for the real MR experiments. Fully sampled k-space data, its inverse DFT as a ground truth, the undersampling mask, and the undersampled k-space data.}\label{RealDataSet}
\end{figure}

\begin{table}[tp!]
\centering
\caption{Comparison of signal-to-noise ratio, high frequency error norm, iteration numbers, and CPU time for the real MR experiments. MaxIter denotes that the iteration reached the maximum number without satisfying the stopping criterion.}\label{RealTable}
\vspace{-0.20cm}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
\hline
\multirow{2}{*}{Indices}&\multirow{2}{*}{Zero fill}&\multirow{2}{*}{TV \cref{TVModel}}&\multirow{2}{*}{Haar \cref{HaarModel}}&\multirow{2}{*}{TLMRI \cref{TLMRI}}&\multicolumn{3}{|c|}{GIRAF \cref{Schatten}}&\multirow{2}{*}{DDTF \cref{ProposedCSMRIModel}}\\ \cline{6-8}
&&&&&$p=0$&$p=0.5$&$p=1$&\\ \hline\hline
SNR&$11.09$&$14.27$&$14.31$&$\textbf{16.36}$&$15.61$&$15.50$&$13.55$&$16.22$\\ \hline
HFEN&$0.6050$&$0.3477$&$0.3393$&$\textbf{0.2141}$&$0.2788$&$0.2874$&$0.4200$&$0.2456$\\ \hline
NIters&$\cdot$&$277$&$364$&MaxIter&$16$&$18$&$18$&$382$\\ \hline
Time($\mathrm{s}$)&$\cdot$&$3.3$&$22.7$&$50097.5$&$75.3$&$104.3$&$98.3$&$16657.4$\\ \hline
\end{tabular}
\end{table}

\cref{RealTable} summarizes the SNR, the HFEN, the number of iterations and the CPU time of the restoration results. For visual comparisons, the restored images, the error maps, the restored k-space data, and the zoom-in views are presented in \cref{RealResults,RealResultsk,RealResultsErrorMap,RealResultsZoom}, respectively. We can observe that our proposed model \cref{ProposedCSMRIModel} outperforms \cref{TVModel,HaarModel,Schatten} as in the phantom experiment, but \cref{TLMRI} performs slightly better than \cref{ProposedCSMRIModel} in the real MR experiment. Nonetheless, the difference of SNR between \cref{ProposedCSMRIModel,TLMRI} is $0.14\mathrm{dB}$, and visually they appear nearly indistinguishable, as shown in \cref{RealTLMRI,RealDDTF,RealTLMRIZoom,RealDDTFZoom}. In addition, even though our approach is relatively time-consuming (approximately $4\mathrm{hrs}$ and $30\mathrm{mins}$) compared to \cref{TVModel,HaarModel,Schatten}, its CPU time is far shorter than \cref{TLMRI} which spends $14\mathrm{hrs}$ but fails to meet stopping criterion. Hence, we can conclude that our proposed approach \cref{ProposedCSMRIModel} performs at least comparably to the TLMRI model \cref{TLMRI} with shorter computational times than \cref{TLMRI}.

Finally, we observe that \cref{ProposedCSMRIModel} for the real MR experiment requires $382$ iterations to meet the stopping criterion, which is more than the piecewise constant phantom image, as shown in \cref{PhantomTable,RealTable}. It should be noted that the original k-space data originally contains the thermal noise. It is also possible that compressing a multi coil k-space data into a virtual single coil data generates an error in the k-space data. In any case, it might fail to meet the underlying assumption of the proposed model; the fully sampled k-space data to be restored is modeled by \cref{MRForward} with the piecewise constant proton density $u$ in \cref{uModel}. In fact, since the proposed data driven tight frame based approach is a relaxation of a low rank (two-fold) Hankel matrix completion, the future work will require the analysis on the convergence rate of \cref{Alg1}, which is likely to depend on the incoherence bound (e.g. \cite{J.F.Cai2019}) of the two-fold Hankel matrix constructed by the original k-space data.

\begin{figure}[tp!]
\centering
\hspace{-0.1cm}\subfloat[Fully sampled]{\label{RealOriginal2}\includegraphics[width=3.00cm]{RealOriginal.pdf}}\hspace{0.005cm}
\subfloat[Zero fill]{\label{RealZeroPad}\includegraphics[width=3.00cm]{RealZeroPad.pdf}}\hspace{0.005cm}
\subfloat[TV \cref{TVModel}]{\label{RealTV}\includegraphics[width=3.00cm]{RealTV.pdf}}\hspace{0.005cm}
\subfloat[Haar \cref{HaarModel}]{\label{RealHaar}\includegraphics[width=3.00cm]{RealHaar.pdf}}\hspace{0.005cm}
\subfloat[TLMRI \cref{TLMRI}]{\label{RealTLMRI}\includegraphics[width=3.00cm]{RealFRIST.pdf}}\vspace{-0.20cm}\\
\subfloat[GIRAF$0$]{\label{RealGIRAF0}\includegraphics[width=3.00cm]{RealGIRAF0.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$0.5$]{\label{RealGIRAFHalf}\includegraphics[width=3.00cm]{RealGIRAFHalf.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$1$]{\label{RealGIRAF1}\includegraphics[width=3.00cm]{RealGIRAF1.pdf}}\hspace{0.005cm}
\subfloat[DDTF \cref{ProposedCSMRIModel}]{\label{RealDDTF}\includegraphics[width=3.00cm]{RealDDTF.pdf}}\vspace{-0.20cm}
\caption{Visual comparisons of each restoration models for the real MR experiments.}\label{RealResults}
\end{figure}

\begin{figure}[tp!]
\centering
\hspace{-0.1cm}\subfloat[Sample mask]{\label{RealMask2}\includegraphics[width=3.00cm]{RealMask.pdf}}\hspace{0.005cm}
\subfloat[Zero fill]{\label{RealZeroPadError}\includegraphics[width=3.00cm]{RealZeroPadError.pdf}}\hspace{0.005cm}
\subfloat[TV \cref{TVModel}]{\label{RealTVError}\includegraphics[width=3.00cm]{RealTVError.pdf}}\hspace{0.005cm}
\subfloat[Haar \cref{HaarModel}]{\label{RealHaarError}\includegraphics[width=3.00cm]{RealHaarError.pdf}}\hspace{0.005cm}
\subfloat[TLMRI \cref{TLMRI}]{\label{RealTLMRIError}\includegraphics[width=3.00cm]{RealFRISTError.pdf}}\vspace{-0.20cm}\\
\subfloat[GIRAF$0$]{\label{RealGIRAF0Error}\includegraphics[width=3.00cm]{RealGIRAF0Error.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$0.5$]{\label{RealGIRAFHalfError}\includegraphics[width=3.00cm]{RealGIRAFHalfError.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$1$]{\label{RealGIRAF1Error}\includegraphics[width=3.00cm]{RealGIRAF1Error.pdf}}\hspace{0.005cm}
\subfloat[DDTF \cref{ProposedCSMRIModel}]{\label{RealDDTFError}\includegraphics[width=3.00cm]{RealDDTFError.pdf}}\vspace{-0.20cm}
\caption{Comparisons of error maps for the phantom experiments.}\label{RealResultsErrorMap}
\end{figure}

\begin{figure}[tp!]
\centering
\hspace{-0.1cm}\subfloat[Fully sampled]{\label{RealOriginalk2}\includegraphics[width=3.00cm]{RealOriginalk.pdf}}\hspace{0.005cm}
\subfloat[Zero fill]{\label{RealUndersamplek}\includegraphics[width=3.00cm]{RealUndersamplek.pdf}}\hspace{0.005cm}
\subfloat[TV \cref{TVModel}]{\label{RealTVk}\includegraphics[width=3.00cm]{RealTVk.pdf}}\hspace{0.005cm}
\subfloat[Haar \cref{HaarModel}]{\label{RealHaark}\includegraphics[width=3.00cm]{RealHaark.pdf}}\hspace{0.005cm}
\subfloat[TLMRI \cref{TLMRI}]{\label{RealTLMRIk}\includegraphics[width=3.00cm]{RealFRISTk.pdf}}\vspace{-0.20cm}\\
\subfloat[GIRAF$0$]{\label{RealGIRAF0k}\includegraphics[width=3.00cm]{RealGIRAF0k.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$0.5$]{\label{RealGIRAFHalfk}\includegraphics[width=3.00cm]{RealGIRAFHalfk.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$1$]{\label{RealGIRAF1k}\includegraphics[width=3.00cm]{RealGIRAF1k.pdf}}\hspace{0.005cm}
\subfloat[DDTF \cref{ProposedCSMRIModel}]{\label{RealDDTFk}\includegraphics[width=3.00cm]{RealDDTFk.pdf}}\vspace{-0.20cm}
\caption{Comparisons of k-space data for the phantom experiments in the log scale. All restored k-space data are displayed in the window level $[0,9]$ for fair comparisons.}\label{RealResultsk}
\end{figure}

\begin{figure}[tp!]
\centering
\hspace{-0.1cm}\subfloat[Fully sampled]{\label{RealOriginalZoom}\includegraphics[width=3.00cm]{RealOriginalZoom.pdf}}\hspace{0.005cm}
\subfloat[Zero fill]{\label{RealZeroPadZoom}\includegraphics[width=3.00cm]{RealZeroPadZoom.pdf}}\hspace{0.005cm}
\subfloat[TV \cref{TVModel}]{\label{RealTVZoom}\includegraphics[width=3.00cm]{RealTVZoom.pdf}}\hspace{0.005cm}
\subfloat[Haar \cref{HaarModel}]{\label{RealHaarZoom}\includegraphics[width=3.00cm]{RealHaarZoom.pdf}}\hspace{0.005cm}
\subfloat[TLMRI \cref{TLMRI}]{\label{RealTLMRIZoom}\includegraphics[width=3.00cm]{RealFRISTZoom.pdf}}\vspace{-0.20cm}\\
\subfloat[GIRAF$0$]{\label{RealGIRAF0Zoom}\includegraphics[width=3.00cm]{RealGIRAF0Zoom.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$0.5$]{\label{RealGIRAFHalfZoom}\includegraphics[width=3.00cm]{RealGIRAFHalfZoom.pdf}}\hspace{0.005cm}
\subfloat[GIRAF$1$]{\label{RealGIRAF1Zoom}\includegraphics[width=3.00cm]{RealGIRAF1Zoom.pdf}}\hspace{0.005cm}
\subfloat[DDTF \cref{ProposedCSMRIModel}]{\label{RealDDTFZoom}\includegraphics[width=3.00cm]{RealDDTFZoom.pdf}}\vspace{-0.20cm}
\caption{Zoom-in views of \cref{RealResults}. The yellow arrows indicate the regions worth noticing.}\label{RealResultsZoom}
\end{figure}

\section{Conclusion and future directions}\label{Conclusion}

In this paper, we propose a new off-the-grid regularization based CS-MRI reconstruction model for the piecewise constant image restoration in the two dimensional FRI framework \cite{G.Ongie2016}\cite{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 sparsely represent a given image. 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.

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 to the analysis on the data driven tight frame, we also need to analyze the relation between the convergence rate of \cref{Alg1} and the incoherence bound of the underlying original k-space data. 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 \cref{uModel}.

Finally, to broaden the scope of applications, it is also likely to extend the idea in this paper to 1) the piecewise smooth image restoration framework, such as the total generalized variation model \cite{K.Bredies2010} and the combined first and second order TV model \cite{M.Bergounioux2010,K.Papafitsoros2014}, by considering the higher order derivatives; 2) the ``blind'' multi-coil CS-MRI reconstruction task, whose inverse problem corresponds to a multi-channel blind deconvolution problem with the unknown coil sensitivity map. For the blind multi-coil CS-MRI reconstruction, even though it is not clear for us at this moment how to extend our approach, this is nevertheless definitely a future direction we would like to work on. Motivated by the recent works on the deep learning techniques combined with the structured low rank matrix approaches \cite{M.Jacob2019,T.H.Kim2019,A.Pramanik2019a,A.Pramanik2019}, we are also interested in developing a deep learning framework based on our approach.

redpajama_set_name: RedPajamaArXiv


Join the Komodo Mailing List to be first to receive major announcements and exclusive.
Previous Story Haryana enhances ex-gratia for kin of police.
This channel is dedicated for you - our community members - to discuss with each about coin related issues.
Get a constantly updating feed of breaking news, fun stories, pics, memes, and videos just for you. We are making this announcement internally at Coinbase and to the public at the same time to remain transparent with our customers about support for future assets.
The messaging app is counting on using its Kin cryptocurrency for.Getting your bitcoin wallet address. Karnataka, INDIA - Duration:.Check Cryptocurrency market prices at WeTalkCoins. in the wake of the very recent announcement of a listing for Ethereum.
The announcement of the state government to give relief and job.Kin Will Continue to Use Ethereum, While Running on Stellar in Parallel for Digital Service Transactions.
Verge Currency (XVG) is a secure and anonymous cryptocurrency, built with a focus on privacy. Komodo Platform is excited to announce that KMD coins can now be.
March 07, 2018 22:26. RS Coin. Neo. 26.08.2018: FOAM Protocol. Ethereum. 6: 31.07.2018: View all ICOs.Torrential rains in many parts of Karnataka have resulted in many dams,.
Women over 45 years of age may soon go on Haj without male kin. Naqvi made the announcement after a Haj review committee. more from times of india News.