Compressed sensing MRI: a review from signal processing perspective

MR forward models

Before we introduce the concept of compressed sensing MR, we begin with the discussion of MR forward model. Mathematically, the forward model for k-space measurement can be described by.

where r ∈ R^d and k ∈ R^d, d = 2, 3 denote the image domain and k-space coordinates, respectively, C is the number of the coil and the i-th coil image γⁱ is given by

Here, ν(r) denotes the contrast weighted bulk magnetization distribution, and sⁱ(r) is the corresponding coil sensitivity of the i-th coil. Although the expression may give an impression that the measurement is two- or three- dimensional, this is indeed originated from a 1-D measurement since the k-space trajectory is a function of time, i.e. k:= k(t), and we acquire one k-space sample at each time point t. This makes the MR imaging inherently slow, since we need to scan through a 3-D object via 1-D trajectories.

Dynamic MRI is another important MR technique to monitor dynamic processes such as brain hemodynamics and cardiac motion. Among the various forms of dy- namic MR modeling, here we mainly focus on the k-t formulation. Specifically, consider a discrete imaging equation for cartesian trajectory for simplicity. Because the samples along the readout direction are fully sampled, most of the dynamic MR formulation is applied separably after taking the Fourier transform along the read- out direction. More specifically, let γ(s, t) denote the unknown image content (for example, proton density, T1/T2 weighted image, etc.) on the spatial coordinate s along the phase encoding line at time instance t. Then, the k-t space measurement b(k, t) at time t is given by

where γ(s, t) is the spatio-temporal image which may be weighted by coil sensitivity map for the case of parallel imaging.

Throughout the paper, for simplicity we often use the operator notation for (1):

where [Sⁱ] is a diagonal operator comprised of the i-th coil sensitivity. Similarly, we use operator notation for (3).

for the case of parallel imaging.

Compressed sensing theory

Optimization approaches

In practice, the following unconstrained form of optimization problem is often used:

$$ \underset{x}{\min}\frac{1}{2}{\left\Vert \mathrm{y}-\mathrm{Ax}\right\Vert}^2+\uplambda {\left\Vert \Psi \mathrm{x}\right\Vert}_1 $$

(8)

where y is noisy measurement, λ is the regularization parameter, and Ψ refers to an analysis transform such that Ψx becomes sparse. Note that for the special case of Ψ = I, (8) is reduced to (P1).

One of the technical issues in solving (8) is that the cost function is not smooth at Ψx = 0, which makes the corresponding gradient ill-defined. This leads to the development of various techniques, which are culminated as the new type of convex optimizaton theory called proximal optimization [28]. For example, the popular op- timization methods such as forward-backward splitting (FBS) [29], split Bregman iteration [30], alternating directional method of multiplier (ADMM) [31], Douglas- Rachford splitting algorithm [32], and primal-dual algorithm [33] have been devel- oped to solve compressed sensing problems. However, the comprehensive coverages of these techniques deserves another review paper, so in this section we mainly review the ADMM algorithms which have been extensively used for compressed sensing MRI ever since its first introduction to compressed sensing MRI [34].

More specifically, we convert the problem (8) to the following constraint problem:

$$ \underset{\gamma, \mathrm{u}}{\min }\ \frac{1}{2}{\left\Vert y- A\gamma \right\Vert}^2+\uplambda {\left\Vert u\right\Vert}_1 $$

(9)

subject to

$$ \mathrm{u}=\Psi \upgamma . $$

(10)

Then, the associated ADMM is given by

$$ {\gamma}^{\left(k+1\right)}=\mathit{\arg}\ \underset{x}{\min }\ \frac{1}{2}{\left\Vert y- A\gamma \right\Vert}^2+\frac{\upmu}{2}{\left\Vert {\upzeta}^{(k)}+\Psi \upgamma -{\mathrm{u}}^{(k)}\right\Vert}^2 $$

$$ {u}^{\left(k+1\right)}=\mathit{\arg}\ \underset{u}{\min}\kern0.50em \uplambda {\left\Vert u\right\Vert}_1+\frac{\upmu}{2}{\left\Vert {\upzeta}^{(k)}+\Psi {\upgamma}^{\left(k+1\right)}-\mathrm{u}\right\Vert}^2 $$

$$ {\upzeta}^{\left(k+1\right)}={\upzeta}^{(k)}+{x}^{\left(k+1\right)}-{u}^{\left(k+1\right)}, $$

Now, each step of ADMM has a closed-form solution. Algorithm 1 summarises the resulting ADMM iteration, where shrink₁ denotes the soft-tresholding:

$$ {shrink}_1\left(\mathrm{x},\uptau \right)=\operatorname{sgn}\left(\mathrm{x}\right)\ \max \left\{0,|\mathrm{x}|-\uptau\ \right\}, $$

(11)

for the threshold value τ > 0.

In addition to the analysis prior in (8), the total variation (TV) penalty has been also extensively used for imaging applications, because the finite difference operator

can sparsify smooth images. Specifically, the TV minimisation problem assumes the following form:

$$ \underset{x}{\min}\frac{1}{2}{\left\Vert y- Ax\right\Vert}^2+\uplambda \mathrm{TV}\left(\mathrm{x}\right), $$

(12)

where T V (x) is given by

$$ \mathrm{TV}\left(\mathrm{x}\right)={\left\Vert \nabla x\right\Vert}_{1,p}, $$

(13)

where ‖∙‖_{1, p}, p = 1, 2 deotes the l₁/l_p-mixed norm. In d-dimensional space (e.g. d = 2 for images), the discretized implementation can be defined as

$$ {\left\Vert \nabla x\right\Vert}_{1,p}={\sum}_{i=1}^n{\left\Vert \nabla x(i)\right\Vert}_p $$

where ∇x(i) ∈ R^d denotes the gradient of x at the i-th coordinate and n denotes the number of discretizated samples.

In order to apply ADMM for (12), we need to focus on the primal formulation of the total variation penalty:

$$ TV(x)={\left\Vert \nabla x\right\Vert}_{1,p}={\sum}_{i=1}^n{\left\Vert \nabla x(i)\right\Vert}_p $$

(14)

Now, we define a splitting variable u(i) = ∇x(i) ∈ R^d. Then, the constraint opti- mization formulation is given by

$$ \underset{x,{\left\{u(i)\right\}}_{i=1}^n}{\min}\frac{1}{2}{\left\Vert y- Ax\right\Vert}^2+\uplambda {\sum}_{i=1}^n{\left\Vert u(i)\right\Vert}_p $$

(15)

$$ \mathrm{Subject}\ \mathrm{to}\ \mathrm{u}\left(\mathrm{i}\right)=\nabla x(i),i=1,\cdots, n $$

(16)

Then, the associated ADMM is given by

$$ {x}^{\left(k+1\right)}=\mathit{\arg}\ \underset{x}{\min }\ \frac{1}{2}{\left|\left|y- Ax\right|\right|}^2+\frac{\upeta}{2}{\sum}_i{\left\Vert {\upzeta}^{(k)}(i)+\nabla x(i)-{\mathrm{u}}^{(k)}(i)\ \right\Vert}^2 $$

$$ {u}^{\left(k+1\right)}(i)=\mathit{\arg}\ \underset{u(i)}{\min}\kern0.50em \uplambda {\left\Vert u(i)\right\Vert}_p+\frac{\upeta}{2}{\left\Vert {\upzeta}^{(k)}(i)+\nabla {x}^{\left(k+1\right)}(i)-\mathrm{u}\left(\mathrm{i}\right)\ \right\Vert}^2 $$

$$ {\upzeta}^{\left(k+1\right)}(i)={\upzeta}^{(k)}(i)+\nabla {x}^{\left(k+1\right)}(i)-{u}^{\left(k+1\right)}\left(\mathrm{i}\right) $$

Each step has the closed form expression. Algorithm 2 summarises the TV-ADMM algorithm, where shrink_vec,2(x) for x ∈ Rd denotes the vector shrinkage [34]:

$$ {shrink}_{vec,2}\left(x,\uptau \right)=\frac{x}{\left|x\right|}\max \left\{0,\left|x\right|-\uptau\ \right\},x\in {\mathrm{R}}^d $$

(17)

for the threshold value τ > 0.

Basic MR ingredients for compressed sensing

In order to apply compressed sensing theory for specific imaging applications, the unknown signal should be sparse in some transform domain, and the sensing matrix should be sufficiently incoherent. This section shows why these conditions can be readily satisfied in MR imaging, which is one of the main reasons that allows for successful applications of CS theory to MR imaging.

Sparsity of MR images

An MR image is rarely sparse in its own. However, one of the important observations that led to the successful applications of CS to MR is that the sparsity is closely related to signal redundancies. This is because redundant signals can be easily converted to sparse signals using some transforms.

Basically, there are three major directions that have been investigated in com- pressed sensing MRI: 1) the spatial domain redundancy, 2) the temporal domain redundancy, and 3) the coil domain redundancy. For example, as shown in Fig. 1(a), natural images can be sparsely represented in finite difference or wavelet transform domain, although the image is not sparse in its own. This observation is the main idea that allows for total variation (TV) and wavelet transform approaches for image denoising, and reconstruction. Accordingly, TV and wavelets have been the main transforms that have been extensively used in most of the CS MRI researches. On the other hand, dynamic MR images such as cardiac cine, functional MRI, and MR parameter mapping have significant redundancy along the temporal dimension as shown in Fig. 1(b). For example, if the image is perfectly periodic, then temporal Fourier transform may be the optimal transform to sparsify the signal. However, in many dynamic MR problems, the temporal variations are dependent on the MR physics as well as specific motion of organs, so the analytic transform such as Fourier transform may not be an optimal solution, but more data-driven approaches such as PCA or dictionary learning are better options. Indeed, these observation leads to dictionary learning approaches that will be discussed later.

On the other hand, the coil redundancy is somewhat distinct compared to the spatial- and temporal- domain redundancy. As shown in Fig. 1(c), the redundancy of the coil image stems from the underlying common images, which results in cross- channel redundancies:

$$ {s}^j(r){\gamma}^i(r)={s}^i(r){\gamma}^j(r),\mathrm{i}\ne \mathrm{j},\forall \mathrm{r} $$

(18)

where sⁱ denotes the i-th coil sensitivity map and γⁱ is the coil image. The rela- tionship in (18) is easy to show since the coil image is given by (2). The main idea of the classical parallel MRI is to exploit this relationship. Specifically, the sensi- tivity encoding (SENSE) [1] exploits the image domain redundancy described in (18), whereas the k-space domain approaches such as GRAPPA [3] exploits its dual relationship in the k-space:

$$ {\widehat{s}}^j(r)\ast {\widehat{\gamma}}^i(r)={\widehat{s}}^i(r)\ast {\widehat{\gamma}}^j(r),\mathrm{i}\ne \mathrm{j},\forall \mathrm{r} $$

(19)

where ∗ denotes the convolution, and \( {\widehat{s}}^j \), \( {\widehat{\gamma}}^j \)enote the Fourier transform of s^jand γ^jrespectively. Later we will describe how these two expressions of the coil dimensional redundancies have been exploited in compressed sensing MRI.

Incoherent sampling pattern

In compressed sensing MRI, downsampling pattern is very important to impose the incoherence, but its realization is limited by MR physics. For example, in the 2-D acquisition, the readout direction should be fully sampled, so there is only freedom along the phase encoding direction in designing incoherence sampling patterns. Figure 2(a)-(c) shows the examples of realizable sampling patterns that have been exploited in the literature: a) cartesian undersampling, b) radial trajectory, and c) spiral trajectory. In particular, the radial and spiral trajectories, which have been studied even before the advanced of the compressed sensing [35], have more incoherent radial sampling patterns compared to the cartesian undersampling.

On the other hand, if we deal with 3-D imaging or dynamic imaging, there are more rooms for sampling pattern design, since there are two dimensional degree of freedom. Figure 2(d) shows an example of 2-D random sampling pattern. For the case of wavelet transform-based sparsity imposing prior, Lustig et al. [14] showed that the incoherence can be optimized by generalizing the notion of PSF to Transform Point Spread Function (TPSF). Specifically, TPSF measures how a single transform coefficient of the underlying object ends up influencing other transform coefficients of the measured undersampled object. To calculate the TPSF, a single point in the wavelet transform space at the i-th location is transformed to the image space and then to the Fourier space. Once the corresponding Fourier space data is subsampled by the given downsampling pattern, its influence in the wavelet transform domain can be calculated by taking inverse transform followed by the inverse wavelet trans- form. To have the best incoherence property, the side of TPSF should be as small as possible, such that the side lobe contribution can be removed by shrinkage op- eration. This was used as the main criterion for sampling pattern design.

Non-Cartesian compressed sensing MRI

Compressed sensing with non Cartesian sampling has been also extensively stud- ied, since they are really a great combination given the sampling behavior of non Cartesian sampling schemes and the incoherence requirement in compressed sens- ing reconstruction. Aside from the radial and spiral sampling patterns discussed before, the works in [71–83] have fo- cused on designing better sampling trajectories with good incoherent properties. For example, Haldar et al. [82] and Puy et al. [83] proposed random phase encod- ing scheme to maximize the incoherency of the sensing matrix. However, one of the main technical issues associated with non-cartesian compressed sensing MRI is that the fast reconstruction trick shown in (21) cannot be used, which increases the overall computational time.

Combination of parallel imaging with CS

Recall that parallel MRI (pMRI) [1, 3] exploits the diversity in the receiver coil sen- sitivity maps that are multiplied by an unknown image. This provides additional spatial information for the unknown image, resulting in accelerated MR data acquisition through k-space sample reduction. Because the aim of the parallel imaging and CS approaches is similar, extensive research efforts have been made to syner- gistically combine the two for further acceleration [18, 20, 38, 40, 84].

One of the most simplest approaches can be a SENSE type approach that explic- itly utilizes the estimated coil maps to obtain an augmented compressed sensing problem [18, 20, 38, 84]. More specifically, if the coil sensitivity is known and given by the sensitivity[Sⁱ], i = 1, ⋯, C,then the SENSE type compressed sensing MRI problem can be formulated as

The optimization framework is the standard optimization framework under sparsity constraint, so proximal optimization algorithms can be used to solve this problem.

On the other hand, l₁-SPIRiT (l₁-iTerative Self-consistent Parallel Imaging Re- construction) [40] utilizes the GRAPPA type constraint as an additional constraint for a compressed sensing problem:

where M is an image domain GRAPPA operator. In both approaches, an accurate estimation of coil sensitivity maps or GRAPPA kernel is essential to fully exploit the coil sensitivity diversity. To address this problem, Uecker et al. [85] developed a novel eigen space method to extract the coil sensitivity maps directly from the k-space data, which is one of the most popular methods widely used by MR researchers.

Blind compressed sensing MR using dictionary learning

Blind compressed sensing approaches attempted to simultaneously reconstruct the underlying image as well as the sparsifying transform from highly undersampled measurements. Ravishankar and his colleagues pioneered two distinct approaches - synthesis dictionary learning [48] and analysis transform learning [86] - when the underlying sparsifying transform is unknown a priori.

k-t SPARSE

The k-t SPARSE by Lustig et al. [88] is an earliest version of compressed sensing dynamic MRI. Recall that the k-space measurement b(k, t) at time t is given by

where ρ(s, f) denotes the temporal Fourier transform of γ(s, t).

Note that p(s, f) is usually sparse because the periodic motions from heart or slow varying motions from fMRI can be easily sparsified using the temporal Fourier transform. Accordingly, the k-t data can be represented as mapping from the spatial- temporal image:

where A: = DF and D is a k-t downsampling pattern, and F now becomes a 2D Fourier transform. In addition to the temporal Fourier transform, the authors in [88] used the wavelet transform in the spatial dimension to exploit the spatial redundancy. Then, k-t SPARSE is formulated based on the following optimization problem:

where W denotes the spatial wavelet transform, respectively.

K-t FOCUSS

Preliminary CS dynamic MRI approaches [36, 88] were seemingly different from the classical k-t approach such as k-t BLAST/SENSE [5]. One of the most impor- tant contributions of the k-t FOCUSS by Jung et al. [18, 20] was to reveal that the compressed sensing dynamic MRI is not very different from the classical k-t approaches, but rather it can be obtained by a very simple modification of the classical k-t BLAST/SENSE to ensure significant performance improvement.

More specifically, rather than using wavelet transform, k-t FOCUSS exploited the x-f domain sparsity. Then, a standard compressed sensing formulation would be:

to enforce the sparseness in x-f image ρ. However, there exists two additional novel- ties in the k-t FOCUSS. First, rather than directly enforcing the sparseness of the x-f image, the k-t FOCUSS further sparsifies the x-f image using the initial estimate.

Specifically, let ρ₀ be the predictable initial estimate of ρ. Then, the residual

should be much more sparse than the original spatio-temporal image. Second, rather than directly using the l₁ minimization, k-t FOCUSS employed the reweighted norm approaches. This results in the following equivalent minimization problem:

Then, the normal equations with respect to x and v are given by

One of the most powerful observations in [18, 20] was that the first iteration of (38) has very similar form to the classical k-t BLAST/SENSE algorithm, except the power factor of weighting matrix. This observation led to an innovative idea to con- vert k-t BLAST/SENSE to compressed sensing approach. More specifically, by using incoherence sampling patterns, multiple iterations and correct weighting factor for the diagonal matrix, the authors of k-t FOCUSS [18, 20] clearly demonstrated the performance improvement. This observation suggested that the improvement by the classical k-t algorithms such as k-t BLAST/SENSE [5] was not from the Bayesian perspective as the original authors of [5] had claimed, but indeed is originated from exploiting the sparsity int the spatio-temporal domain [18, 20]. Furthermore, by simply modifying the weight factor and sampling patterns, several additional iter- ation can significantly improve the performance of k-t BLAST/SENSE.

Another powerful aspect of k-t FOCUSS was that the idea can be easily extended to exploit the sparsity in other transform domains. For example, the residual step in (36) can be interpreted as sparsity promoting step by subtracting the temporal mean images. Thus, Jung et al. [20] proposes motion estimated and compensated modification of k-t FOCUSS to make the residual signal much sparser. More specifi- cally, rather than subtracting the temporal mean values, they subtracted the motion estimated frame. Note that motion estimation and compensation (ME/MC) is an essential step in video coding that uses motion vectors to exploit the temporal re- dundancies between frames [89, 90]. In order to employ ME/MC within dynamic MRI, there are several technical issues to address. First, at least one reference frame is required. This issue can be easily resolved if we acquire fully sampled data in one frame as often done in dynamic MRI. The main technical difficulty, however, comes from the existence of the low quality current frame. Fortunately, this issue can be addressed using an additional reconstruction step before the ME/MC [20].

where D denotes the learned temporal dictionary based from the images, whereas c denotes the coefficients. Then, the imaging problem can be formulated as

For example, in order to find the temporal basis that can sparsify ρ, Jung et al. [18, 20] performed the the singular value decomposition (SVD) after the image reconstruction, then the new dictionary is used to estimate the new coefficients. This procedure is closely related to the partial separable function (PSF) [91] and k-t SLR (k-t sparse and low-rank decomposition) [49], which will be reviewed soon.

Partially separable function (PSF) approach

In the PSF model by Liang et al. [91], the k-t samples are assumed to be decomposed in the following form

has at most rank L.

In dynamic CS MRI, many of the k-t samples are missing and we are interested in finding the missing components. Hence, by utilizing the low-rankness of the Ca-sorati matrix, the missing k-t samples can be estimated using a low rank matrix completion algorithm. In particular, the authors in [92, 93] proposed the following matrix factorization approach:

where A denotes the k-t sampling operator that indicates the missing samples by 0, and \( U\in {\mathrm{C}}^{m\times \mathrm{L}} \), \( V\in {\mathrm{C}}^{m\times \mathrm{L}} \) is used for the low rank matrix factorization B = U V ^H, and the optimization is performed for U and V alternatingly by fixing the other matrix using the previous estimate. Because the exact rank of B is not known, the authors proposed the incremented powerFactorization (IRFP) algorithm [92], where (40) starts with L = 1 by increasing order with the initialization of the power factorization of L + 1 from that of L. To avoid an overfitting, the algorithm is terminated as soon as the data fidelity is below some threshold values.

K-t SLR: K-t sparse and low rank approach

k-t Sparse and Low Rank Approach (k-t SLR) model by Lingala et al. [49] is a more systematic way of learning both basis and sparse coefficients. In this approach, the spatio-temporal signal ρ(x, t) is first rearranged in a matrix form

Then, using the low rank prior, the optimisation problem can be formulated as

\( \underset{\Gamma}{\min }{\left\Vert b-A\left(\Gamma \right)\right\Vert}^2+\uplambda \upvarphi \left(\Gamma \right) \),

where A: = DF is now a downsampled Fourier transform. Here, the rank prior is approximated using the general class of Schatten p-functionals, specified by

where {σ_i} denotes the singular values of Γ.

Basic theory

Compared to the standard compressed sensing approaches, k-space structured low- rank approaches such as SAKE [55], LORAKS [56], ALOHA [57, 59, 60] and GIRAF [61] are relatively new, but has significant potentials in MRI imaging. These ap- proaches are all derived by the k-space convolution relationship and can be used for both static and dynamic imaging. So we first discuss a matrix representation of the convolution. For simplicity, we will consider the 1-D notation.

where γⁱ denotes the unknown i-th coil images and yⁱ corresponds to its k-space data. Note that the matrix representation of a k-space convolution of yⁱ with a d-tap filter h is given by

where \( \mathcal{C}\left({y}^i\right) \) denotes the convolution matrix constructed by the vector yⁱ:

By defining Y = [y¹, ⋯, y^C], we can further defined the extended Hankel matrix

In the following, we will explain how these Hankel structured matrix has been utilised for accelerated MRI.

SAKE

A calibrationless parallel imaging reconstruction method, termed simultaneous au- tocalibrating and k-space estimation (SAKE), is a data-driven, coil-by-coil recon- struction method that does not require a separate calibration step for estimating coil sensitivity information [55]. SAKE is based on the following observation in GRAPPA:

which implies that the i-th k-space measurement can be represented as the linear combination of the filtered k-space data from other coils. In matrix form, this recursive relationship implies the existence of the null space of the Hankel matrix.

In other word, H_d|C(Y) is low-ranked. Therefore, SAKE solves the following low- rank matrix completion problem to interpolate the missing k-space data:

where P_Ω(·) denotes the projection on the measured k-space samples on the index set Ω. The problem was solved using the iterative singular value shrinkage method [94].

LORAKS

The low-rank modeling of local-space neighborhoods (LORAKS) [56] was inspired by the finite support condition. More specifically, if the object γ has finite support, we can easily find the function w such that

This results in a convolution relationship in k-space

where y and h denote the Fourier spectrum of γ and w, respectively. Therefore, this gives a single channel version of the low-rank condition, which results in the following rank minimization problem:

ALOHA and GIRAF

Annihilating filter-based low rank Hankel matrix (ALOHA) approach [57–60, 62] and GIRAF (Generic Iterative Reweighted Annihilating Filter) [61] can be considered as the full generalization of SAKE and LORAKS for general class of signals with the finite rate of innovations (FRI) for MR measurements. Moreover, the approaches have unified the parallel imaging and compressed sensing as a k- space interpolation with performance guarantees [62], and can be used for artifact correction [58]. This section describes the fundamental dual relationship between transform domain sparsity and low rankness in reciprocal domain, which is the key ingredient. For better readability, we provide here a high level description by assuming 1-D signals.

The Fourier CS problem of our interest is to recover the unknown signal x(t) from the Fourier measurement:

when the support of the time domain signal x(t) is τ. Then, discrete Fourier data at the Nyquist rate is defined by:

We also define a length (r + 1)-annihilating filter hˆ[k] for xˆ[k] that satisfies

The existence of the minimum length finite length annihilating filter has been ex- tensively studied for FRI signals [95–97]. Let r + 1 denotes the minimum size of annihilating filters that annihilates discrete Fourier data \( \widehat{\mathrm{x}}\left[\mathrm{k}\right] \). Then, a d-tap anni- hilating filter h with d > r + 1 can be easily obtained by convolving an appropriate size FIR filter with the minimum length annihilating filter. In matrix form, this is equivalent to

where the Hankel structure matrix \( {\mathcal{H}}_{\mathrm{d}}\left(\hat{x}\right) \) is constructed as

Assume that min(n − d + 1, d) > r. Then, we can show the following low rank property [62]:

Thanks to the low-rankness of the associated Hankel matrix, the missing k-space data can be easily interpolated using the following low-rank matrix completion [62]:

where P_Ω is the projection operator on the sampling location Ω. Moreover, as shown in [62], the low-rank matrix completion approach (55) does not compromise any optimality compared to the standard Fourier CS.

Note that signals may not be sparse in the image domain, but can be sparsified in a transform domain. In fact, this was the main idea of the compressed sensing. Specifically, the signal x of our interest is a non-uniform spline that can be represented by:

where L denotes a constant coefficient linear differential equation that is often called the continuous domain whitening operator in [98, 99]:

and w is a continuous sparse innovation:

For example, if the underlying signal is piecewise constant, we can set L as the first differentiation. In this case, x corresponds to the total variation signal model. Then,

by taking the Fourier transform of (56), we have

where

Accordingly, the Hankel matrix \( \mathcal{H}\left(\hat{z}\right) \) from the weighted spectrum zˆ(ω) satisfies the following rank condition:

Thanks to the low-rankness, the missing Fourier data can be interpolated using the following matrix completion problem:

or, for noisy Fourier measurements \( \widehat{\mathrm{y}} \),

where ⊙ denotes the Hadamard product, \( \widehat{\mathrm{l}} \) and \( \widehat{\mathrm{x}} \) denotes the vectors composed of full samples of \( \widehat{\mathrm{l}} \)(ω) and \( \widehat{\mathrm{x}} \)(ω), respectively. After solving (P_w), the missing spec- tral data xˆ(ω) can be obtained by dividing by the weight, i.e. \( \widehat{x}\left(\omega \right)=\mathrm{m}\left(\upomega \right)/\widehat{l}\left(\omega \right) \) assuming that \( \widehat{l}\left(\omega \right)\ne 0 \).

The idea can be generalized for any transform domain sparse signals as long as the transform can be represented using shift-invariant filters. Wavelet domain sparse signal belongs to this class. In this case, the weight kernel in the Fourier domain is obtained as the spectrum of the subband filters [57–60, 62]. For example, Fig. 3 showed the construction of the Hankel matrix for the MR parameter mapping, where the k-space weighting from wavelet weighting is applied only along the phase encoding direction, whereas no weighing is applied along the t-domain since the correspond temporal spectrum is already sparse [59].

The idea can be also easily generalised to the parallel imaging by exploiting (19). Specifically, (19) can be equivalently represented using the matrix representation:

This implies that an exten^ded Hankel matrix \( {\mathcal{H}}_{\mathrm{d}\mid \mathrm{C}}\left(\left[{\widehat{\gamma}}^1\cdots {\widehat{\gamma}}^C\right]\right) \) in (46) is low ranked.

For example, the multi-channel construction of Hankel matrix for MR parameter mapping [59] is shown in Fig. 4.

One of the most important advantages of the Hankel matrix formulation is that the coil diversity can be readily exploited in addition to the image domain redundancy in a unified framework. This makes the separate coil sensitivity estimation unnecessary.