RIP Sensing Matrices Construction for Sparsifying Dictionaries with Application to MRI Imaging

Abstract

Practical applications of compressed sensing often restrict the choice of its two main ingredients. They may (i) prescribe the use of particular redundant dictionaries for certain classes of signals to become sparsely represented or (ii) dictate specific measurement mechanisms which exploit certain physical principles. On the problem of RIP measurement matrix design in compressed sensing with redundant dictionaries, we give a simple construction to derive sensing matrices whose compositions with a prescribed dictionary have with high probability the RIP in the $klog(n/k)$ regime. Our construction thus provides recovery guarantees usually only attainable for sensing matrices from random ensembles with sparsifying orthonormal bases. Moreover, we use the dictionary factorization idea that our construction rests on in the application of magnetic resonance imaging, in which also the sensing matrix is prescribed by quantum mechanical principles. We propose a recovery algorithm based on transforming the acquired measurements such that the compressed sensing theory for RIP embeddings can be utilized to recover wavelet coefficients of the target image, and show its performance on examples from the fastMRI dataset.

1. Introduction

Compressed sensing provides a framework under which sparse or compressible signals can be stably reconstructed from far fewer linear measurements than their ambient dimension [1,2]. The number of required measurements depends on the signal complexity in terms of sparsity, and properties of the sensing matrix with respect to sparse vectors, such as the restricted isometry property (RIP) [3]. A sensing matrix S, which embeds high-dimensional signals into a lower-dimensional measurement space, is said to have the RIP of order k if there exists a constant ${\delta }_k\in [0,1)$, such that for any k-sparse vector $x\in {\mathbb{R}}^n$:

$$\begin{array}{c}\hfill (1-{\delta }_k){\parallel x\parallel }_2^2\le {\parallel Sx\parallel }_2^2\le (1+{\delta }_k){\parallel x\parallel }_2^2.\end{array}$$

(1)

Improving upon [4], Candès [5] showed that k-sparse and compressible vectors can be stably recovered from observations $y=Sx+\eta$ with measurement error $\eta$ bounded by $ϵ$, given that S has the RIP with ${\delta }_{2k}<\sqrt{2}-1$ (which has been improved to 0.4652 [6]) via the sparsity promoting convex program:

$$\begin{array}{c}\hfill \underset{x}{minimize}{\parallel x\parallel }_1\mathrm{subject}\mathrm{to}{\parallel y-Sx\parallel }_2\le ϵ.\end{array}$$

(2)

For certain $\lambda >0$, (2) can be equivalently reformulated to the problem of minimizing over x the unconstrained objective $\frac{1}{2}{\parallel y-Sx\parallel }_2^2+\lambda {\parallel x\parallel }_1$ consisting of a data-fidelity and regularization term, which can be solved using iterative methods [7].

While it is NP-hard to verify the RIP for a given matrix, random matrices from sub-Gaussian or Bernoulli ensembles do possess the RIP with high probability—also when composed with orthonormal bases [8]. In practice, however, many classes of signals are sparse only with respect to a redundant dictionary or tight frame that is non-orthonormal (e.g., Gabor, curvelet, wavelet, or data-driven learned dictionaries). In this case, the sensing matrix S in (1) has to be replaced by the composition $SD$, which no longer possesses the RIP. This prevents the RIP as recovery guaranteeing tool in many CS applications.

Contributions: We first present a method to derive a sensing matrix $S\in {\mathbb{R}}^{m\times l}$ for a given sparsifying dictionary $D\in {\mathbb{R}}^{l\times n}$, such that with high probability, $SD$ has the RIP for m on the order of $klog(n/k)$, where k is the sparsity level of a coefficient vector. We use a sufficient condition for the RIP of random matrices that satisfy a concentration of measure inequality [8]. Starting from any random matrix A of the dictionary dimension for which a random row selection $\mathcal{E}A\in {\mathbb{R}}^{m\times n}$ satisfies this concentration inequality, a tailored sensing matrix can be derived from any factorization $D=GAH$, in which $G\in {\mathbb{R}}^{l\times l}$ is invertible and $H\in {\mathbb{R}}^{n\times n}$ is orthonormal: for the sensing matrix $S:=\mathcal{E}{G}^{-1}$, the composition $SD$ has then the desired RIP with high probability. We show in Section 3 that the required factorization exists whenever D and A have equal rank and detail constructions. A particular implication is that one can, thus, obtain the same RIP-based recovery guarantees for sensing matrix compositions with general over-complete dictionaries as for Gaussian or Bernoulli ensembles with orthonormal bases.

We utilize dictionary factorization in the context of CS MRI applications to showcase the viability of our method for a significant purpose. The primary aim of accelerating MRI scans is to lessen the time patients spend in the machines, thereby reducing the discomfort that can result from remaining motionless. The speed at which MRI images are accelerated is affected by the number of rows in an over-complete dictionary. As a result, recovery algorithms aim to use fewer rows while ensuring high-quality image reconstruction. In CS MRI, sensing mechanisms exploit quantum mechanical principles and existing hardware constraints restrict sensing matrix design beyond non-uniformly subsampling the measured (complex-valued) Fourier spatial frequency coefficients of the target image $\tilde{x}=Dx$, which we suppose is synthesized from (real-valued) sparse coefficients x with respect to a dictionary D. The sensing matrix $S=\mathcal{R}F$ is thus modeled as the product of a discrete Fourier transform F and a row subsampling $\mathcal{R}$. The CS recovery of x can then be formulated as the ${\ell }_1$-synthesis problem:

$$\begin{array}{c}\hfill \underset{x}{minimize}{\parallel x\parallel }_1\mathrm{subject}\mathrm{to}{\parallel y-\mathcal{R}FDx\parallel }_2\le ϵ\end{array}$$

(3)

in which $\mathcal{R}F$ possesses the RIP  [9,10,11], while, in general, $\mathcal{R}FD$ does not, such that the RIP recovery guarantees do not directly apply. Using the dictionary factorization idea, we attempt to neutralize the effect of F in order to obtain an RIP synthesis problem. We do so by choosing (real) factors G such that (both real and imaginary parts of) $\mathcal{R}FG$ optimally match $\mathcal{R}$, so that the factorizations $GAH$ that optimally approximate D allow, by virtue of $\mathcal{R}FD\sim \mathcal{R}FGAH\sim \mathcal{R}AH$, to replace the constraint in (3) by ${\parallel y-\mathcal{R}AHx\parallel }_2\le ϵ$. Both the real and the imaginary parts of $\mathcal{R}AH$ are incoherent sampling matrices for sparse signals that possess the RIP. We provide a review of CS application in Section 2 and detail the application of CS factorization method to MRI in Section 4, and report in Section 5 numerical experiments comparing our approach to total-variation based CS MRI [12].

Related work: If the sparsifying dictionary in the general ${\ell }_1$-synthesis approach is a tight frame, i.e., $x=D^{\top }\tilde{x}$, recovery guarantees can be derived by considering the ${\ell }_1$-analysis approach of minimizing over $\tilde{x}$ the objective $\parallel D^{\top }\tilde{x}{\parallel }_1$ subject to the constraint $\parallel y-S\tilde{x}{\parallel }_2\le ϵ$. Recovery guarantees for the ${\ell }_1$-analysis approach are connected to the notion of D-RIP [13], which requires the sensing matrix to satisfy the RIP inequality for all images of k-sparse vectors under a tight frame D. Any RIP matrix satisfies the D-RIP when multiplied by a random sign matrix [14]. Unless D is orthonormal, the geometric structures, properties and empirical performances of the ${\ell }_1$-analysis and ${\ell }_1$-synthesis approaches in general differ [15].

A weaker condition than the RIP that can facilitate recovery guarantees for $SD$ with general dictionary D is the mutual incoherence between the sensing matrix and the dictionary [2,16]. However, sparsity ranges for signals to be recovered, as well as necessary number of measurements are more restrictive than for RIP guarantees, and overcomplete dictionaries with highly coherent columns in general lead to large coherence of their product with sensing matrices. Finally, the nullspace property of a sensing matrix gives a necessary and sufficient condition for stable recovery of sparse signals using the convex optimization [17].

2. Compressed Sensing Applications

Compressed Sensing (CS) has found numerous applications across various fields. One significant area is addressing the future data storage challenges [18]. For example, advancements in communication and signal processing are producing vast amounts of data; 5G networks can reach peak data rates of about 300 Mb per second, while an image of a black hole needs several petabytes for storage. In medical signal processing, the volume of required data is also substantial. CS can enhance wireless communication capabilities in cognitive radio networks and enable spatial compression and projection in wireless sensor networks [19]. Additionally, modern imaging techniques rely on arrays of millions of detector pixels. In contrast, CS is foundational to single-pixel cameras, which utilize mask patterns to filter scenes and capture transmitted intensity measurements with a one-pixel detector. This technology is well-suited for imaging beyond the visible spectrum. It allows for precise timing or depth resolution, making it ideal for applications such as infrared imaging and 3D situational awareness in autonomous vehicles [20].

An important CS application, in which the sensing mechanisms is dictated by hardware constraints, is magnetic resonance imaging (MRI). In MRI a scanner collects Fourier domain measurements of a target medical image. Accelerating the speed of MRI data acquisition by reducing the number of required measurements remains of great interest to the medical community. Inference of the true underlying spatial image can be achieved via several CS strategies when non-uniformly under-sampling below the requirements of the Shannon–Nyquist theory within the maximum frequency spectrum, by incorporating the a priori knowledge about sparsity of medical images in a dictionary transform domain or of its spatial gradients [21]. While deep learning approaches have been a recent center of attention to recover data under-sampled in this way [22], recovery guarantees and interpretability set the CS approach apart from such machine learning techniques [23]. Though state-of-the-art machine learning models for MRI have been validated for clinical interchangeability in four-fold data acquisition acceleration [24], they rely on empirical studies and present challenges concerning, for instance, reconstruction hallucinations that can be problematic for interpreting radiologists [25], unknown training data biases [26], or robustness to distribution shifts between training and application data (e.g., scanning technology [27], target anatomy [28], or acceleration factors [29]). Significant issues concerning the stability of distribution shifts are highlighted in [30]. The authors stress the need to apply a method across different datasets featuring anatomical regions and images collected through different scanning protocols. This approach is crucial for showcasing the method’s versatility and evaluating the model’s resilience to distribution shifts. In contrast, the sparsity-driven CS approach can be fine-tuned to perform close to deep learning methods for MRI acceleration, via unrolling algorithms which consist of a small fraction of the parameters employed by deep learning approaches [23].

3. Sensing Matrix Construction

Signal complexity, in terms of sparsity, determines the amount of possible under-sampling in CS. Constructions utilizing randomness can produce RIP matrices for which the number of required measurements m scales linearly with the sparsity level k of the vector to be recovered. Such matrices can be derived from distributions for which the following concentration of measure inequality, resembling (1), holds, such as for subgaussian or Bernoulli ensembles.

Theorem 1

([8]). Let $0<{\delta }_k<1$ and $A\in {\mathbb{R}}^{m\times n}$ be an iid random matrix. If ${\mathbb{E}\parallel Ax\parallel }_2^2={\parallel x\parallel }_2^2$ for all $x\in {\mathbb{R}}^n$, and for any $ϵ\in (0,1)$ the concentration inequality ${\mathbb{P}(\parallel Ax\parallel }_2^2-{\parallel x\parallel }_2^2{|\ge ϵ\parallel x\parallel }_2^2)\le 2e^{-lc\left(ϵ\right)}$ holds for some $c\left(ϵ\right)>0$ and all $x\in {\mathbb{R}}^n$, then there exists $c_1,c_2>0$ (depending only on ${\delta }_k$) such that, whenever $k\le {\displaystyle \frac{c_{1}m}{log(n/k)}}$, the RIP (1) holds for A with probability at least $1-2e^{-c_{2}m}$.

Random matrices satisfying the concentration inequality are universal with respect to orthonormal bases [8], i.e., the same RIP conclusions hold for their products with unitary matrices. We use this universality to derive sensing matrices for any sparsity-inducing dictionary by random row selection of an invertible transform adapted to the dictionary.

Theorem 2.

Let $D\in {\mathbb{R}}^{l\times n}$ ($l\le n$) be a dictionary, $A\in {\mathbb{E}}^{l\times n}$ and $\mathcal{E}\in {\mathbb{E}}^{m\times l}$ matrices such that $\mathcal{E}A$ satisfies the assumptions of Theorem 1. Suppose the dictionary allows a factorization $D=GAH$ for some invertible G and orthonormal H. Then $S:=\mathcal{E}{G}^{-1}\in {\mathbb{R}}^{m\times l}$ is a sensing matrix for D such that, with probability as in Theorem 1, $SD$ has the RIP whenever $m\gtrsim klog(n/k)$.

In this construction $SD$ equals $\mathcal{E}AH$, which implies the claim due to the above mentioned universality. The existence of a factorization $D=GAH$ as assumed in Theorem 2 requires A and D to have equal rank. We next show that the latter also suffices. Given a full-rank sparsifying dictionary, the sensing matrix construction of Theorem 2 can therefore (with probability one) be carried out, starting from a Gaussian matrix A of the same dimensions as the dictionary, and a random row selection $\mathcal{E}$.

Proposition 1.

Let $A,D\in {\mathbb{R}}^{l\times n}$ with $l\le n$. Then, the following statements are equivalent:

(i) 

There exists an invertible $G\in {\mathbb{R}}^{l\times l}$ and orthonormal $H\in {\mathbb{R}}^{n\times n}$ such that $GD=AH$;

(ii) 

There exists an invertible $G\in {\mathbb{R}}^{l\times l}$ with $GD{D}^{\top }{G}^{\top }=A{A}^{\top }$;

(iii) 

A and D have equal rank.

Proof. 

We show that (iii) implies (ii), and in turn (i). Assuming (iii), $A{A}^{\top }$ and $D{D}^{\top }$ have equal rank. Their respective spectral decompositions $Q_A{\Sigma }_A{Q}_A^{\top }$ and $Q_D{\Sigma }_D{Q}_D^{\top }$ can, thus, be assumed to have zeros in the same positions of the diagonal matrices ${\Sigma }_A$, ${\Sigma }_D$. Replacing those zero diagonal entries by ones to result in ${\Sigma }_A^{\prime }$, ${\Sigma }_D^{\prime }$, the matrix $G:=Q_A{\left({\Sigma }_A^{\prime }{\Sigma }_D^{\prime -1}\right)}^{1/2}{Q}_D^{\top }$ is invertible and:

$$\begin{array}{cc}\hfill GD{D}^{\top }{G}^{\top }& =G{Q}_D{\Sigma }_D{Q}_D^{\top }{G}^{\top }\hfill \\ \hfill & =Q_A{\left({\Sigma }_A^{\prime }{\Sigma }_D^{\prime -1}\right)}^{1/2}{\Sigma }_D{\left({\Sigma }_A^{\prime }{\Sigma }_D^{\prime -1}\right)}^{1/2}{Q}_A^{\top }=A{A}^{\top }.\hfill \end{array}$$

Defining $H:=A^{+}GD+N_A{N}_D^{\top }$, with $A^{+}$ denoting the pseudo-inverse of A, and $N_A$, respectively $N_D$, having pairwise orthonormal columns spanning the nullspaces of A, respectively D, implies:

$$AH=A{A}^{+}GD+A{N}_A{N}_D^{\top }=A{A}^{+}GD=GD$$

since $A{A}^{+}$ is the identity on the range of A, which contains the range of G. Moreover:

$$\begin{array}{cc}\hfill H{H}^{\top }& =A^{+}GD{D}^{\top }{G}^{\top }{\left(A^{+}\right)}^{\top }+N_A{N}_A^{\top }\hfill \\ \hfill & =A^{+}A{A}^{\top }{\left(A^{+}\right)}^{\top }+N_A{N}_A^{\top }=A^{+}A+N_A{N}_A^{\top },\hfill \end{array}$$

which is the identity. Finally, note that (i) implies (iii).    □

An alternative factorization can be derived from orthonormal bases for the ranges and nullspaces of prescribed equal rank $A,D\in {\mathbb{R}}^{l\times n}$ as follows. Let $\left[U_A{N}_A\right],\left[U_D{N}_D\right]$ be orthonormal matrices, where the columns of $U_A,U_D$ span the ranges, and the columns of $N_A,N_D$ the nullspaces, of $A,D$. Extend $D{U}_D,A{U}_A$ to invertible square matrices $\widetilde{D{U}_D},\widetilde{A{U}_A}$, appending columns. Then:

$$\begin{array}{c}\hfill G:=\widetilde{D{U}_D}{\widetilde{A{U}_A}}^{-1}\mathrm{and}H:=U_A{U}_D^{\top }+N_A{N}_D^{\top }\end{array}$$

(4)

are invertible, respectively, orthonormal, and $GAH-D=(GA{U}_A-D{U}_D)U_D^{\top }$ is zero since $GA{U}_A=D{U}_D$.

We next detail a factorization for the practically important case of a prescribed tight frame dictionary D. We may then assume $D{D}^{\top }$ to be the identity.

Theorem 3.

Suppose $A,D\in {\mathbb{R}}^{l\times n}$ have full rank $l\le n$, and D is a tight frame. Then $D=GAH$ with invertible, respectively orthonormal:

$$\begin{array}{c}\hfill G:=O{\left(A{A}^{\top }\right)}^{-1/2}\mathit{and}H:=A^{\top }{G}^{\top }D+N_A{N}_D^{\top },\end{array}$$

where $O\in {\mathbb{R}}^{l\times l}$ is any orthonormal matrix, and the columns of $N_A,N_D\in {\mathbb{R}}^{n\times n-l}$ are any orthonormal bases of the nullspaces of A, respectively D. In particular, if A is also a tight frame, then G is orthonormal.

Proof. 

Since $GA{A}^{\top }{G}^{\top }$ is the identity and $A{N}_A$ the zero operator, a direct calculation shows $GAH=D$. While G is by definition invertible, $D{D}^{\top }$ being the identity implies that $H{H}^{\top }=A^{\top }{\left(A{A}^{\top }\right)}^{-1}A+N_A{N}_A^{\top }$ is the identity.    □

4. Application to Accelerated MRI Imaging

The goal of fast MRI is the recovery of an image $Z\in {\mathbb{R}}^{n_1\times n_2}$ from under-sampled Fourier (k-space) measurements. The discrete Fourier transform of Z is the complex matrix:

$$\begin{array}{c}\hfill \widehat{Z}:=F_{1}Z{F}_2\in {\mathbb{C}}^{n_1\times n_2},\end{array}$$

where $F_1\in {\mathbb{C}}^{n_1\times n_1}$, $F_2\in {\mathbb{C}}^{n_2\times n_2}$ are (symmetric and unitary) discrete Fourier transform matrices. Data acquired via accelerated MRI can be modeled as:

$$\begin{array}{c}\hfill Y:=\mathcal{R}\widehat{Z}\in {\mathbb{C}}^{n_1\times n_2},\end{array}$$

(5)

where $\mathcal{R}$ is a $\{0,1\}$-entry $n_1\times n_1$ diagonal matrix that selects rows of Y corresponding to some under-sampling scheme [31]. The image can be assumed to have a sparse representation:

$$\begin{array}{c}\hfill Z=D_{1}X{D}_2^{\top }\end{array}$$

(6)

in the transform domain of wavelet dictionaries $D_i\in {\mathbb{R}}^{n_i\times N_i}$ ($N_i\ge n_i$, $i=1,2$) that implement a separable bivariate wavelet transform. The coefficient matrix $X\in {\mathbb{R}}^{N_1\times N_2}$ then consists of low-pass channel approximation coefficients, and sparse high-pass channel coefficients that capture edges and singular points in the target image. In this formulation, the CS-MRI problem is then the inverse problem of approximating sparse coefficients X that solve:

$$\begin{array}{c}\hfill Y=\mathcal{R}{F}_1{D}_{1}X{D}_2^{\top }{F}_2,\end{array}$$

(7)

that is, to reconstruct Z (via (6)) from the measurements Y which are in general corrupted by noise.

Being tied to the measurements (7) by physical principles, we first transform the acquired data. Starting from a sensing factorization $D_2=G_{2}A{H}_2$ according to Proposition 1, we consider the transformed observation:

$$\tilde{Y}:=Y{F}_2^{\ast }{\left(G_2^{\top }\right)}^{-1}=\mathcal{R}{F}_1{D}_{1}X{\left(A_2{H}_2\right)}^{\top }.$$

We next construct for $D_1$ two (real) factorizations that are adapted to the real and imaginary parts of $F_1$. First, we let:

$$\begin{array}{c}\hfill G_R:=\underset{G\in {\mathbb{R}}^{n_1\times n_1}\mathrm{invertible}}{argmin}{∥\mathcal{R}(Re\left(F_1\right)G-I)∥}_F^2\end{array}$$

(8)

and, second, we let:

$$\begin{array}{c}\hfill H_R:=\underset{H\in {\mathbb{R}}^{n_2\times n_2}\mathrm{unitary}}{argmin}{∥D_1-G_{R}AH∥}_F^2.\end{array}$$

(9)

Analogously, we define $G_I$ and $H_I$ via solving (8) for the imaginary part of $F_1$, before solving (9) for $G_I$. We then let $H:=H_R+i{H}_I$. For details on numerical solutions of (8) and (9) via the augmented Lagrangian approach we refer to the Appendix A.

Our motivation for desiring in (8) small differences between $\mathcal{R}$ and the real (respectively imaginary) parts of $\mathcal{R}{F}_1{G}_R$ (respectively $\mathcal{R}{F}_1{G}_I$) is that then (9) implies that the following difference is small:

$$\begin{array}{c}\hfill \Phi \left(X\right):=\frac{1}{2}{∥\tilde{Y}-\mathcal{R}AHX{\left(A{H}_2\right)}^{\top }∥}_F^2,\end{array}$$

such that we can attempt to approximate sparse coefficients X by solving the ${\ell }_1$-regularized problem:

$$\begin{array}{c}\hfill \underset{X\in {\mathbb{R}}^{N_1\times N_2}}{minimize}\Phi \left(X\right)+\gamma {\parallel X^H\parallel }_1,\end{array}$$

(10)

where $X^H$ denotes the highpass coefficients of X and $\gamma >0$.

5. Experimental Results

5.1. Sensing Matrix of Sparsifying Dictionaries

Experiments were conducted to assess the sensing matrix derived based on the dictionary factorizations proposed in Section 3 for CDF 9/7 wavelets and dictionaries derived via K-SVD, for which synthetically generated univariate sparse signals are then shown to be recoverable with the guarantees of Theorem 2.

The entries of $A\in {\mathbb{R}}^{l\times n}$ are i.i.d. Bernoulli random variables (where $\pm {\displaystyle \frac{1}{\sqrt{n}}}$ with equal probability) and Gaussian random numbers (where $\mathcal{N}(0,n^{-1})$). D and A have equal ranks. Experiments were performed on sparsifying dictionaries of ${\mathbb{R}}^{128\times 1024}$, including $D_{wavelet}$ and $D_{KSVD}$, which were derived using the K-SVD algorithm, based on the set of training vectors obtained from the $256\times 256$ gray-scale test image “Boat” (parameters of the K-SVD process were set at $K=1024$, $L=64$, $numIteration=50$, $InitializationMethod{=}^{\prime }GivenMatri{x}^{\prime }$, and $errorGoal={10}^{-8}$). The image was divided into overlapping patches of $16\times 16$ pixels with stride 2 (in each dimension) and overlaps of 4 pixels (in each dimension), which resulted in 3721 training patches. As the stride is 2 in each dimension, each patch formed a vector of 128 pixels. After the mean of each patch was normalized to zero, the resulting patches were transformed into vectors (via the vec operation) to form the set of training vectors. The dictionary $D_{wavelet}$ was derived from the CDF 9-7 wavelet. The first 640 columns (i.e., each level has 128 columns and there are 5 levels) in $D_{wavelet}$ were obtained from the first 5 levels (the support of a wavelet at level 6 is larger than the size of an image patch), whereas the remainder were generated randomly. The column norms of the sparsifying dictionaries were normalized. Figure 1 and Figure 2 present the sparsifying dictionaries and corresponding matrices G, derived for various A in accordance with (4). The subgraph associated with G in Figure 2 corresponds to the sensing matrix of D. It is worth noting that the lack of structure in G does not contradict some research in compressed sensing literature that focuses on designing structured sensing matrices. The structured sensing matrix approach aims to efficiently manage an extensive system’s computational resources, such as storage and computational efficiency. In contrast, our research prioritizes creating a sensing matrix that complements a given dictionary, with fewer rows than columns, to maintain the RIP. Exploring the possibility of constructing a structured sensing matrix while retaining the RIP for a given dictionary is an interesting question for future investigation.

Under Theorem 2, the following two optimizations derive the same solution, if $\mathcal{E}{G}^{-1}D=\mathcal{A}H$:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle \underset{x}{min}{\parallel z-\mathcal{E}{G}^{-1}Dx\parallel }_2^2}\hfill \\ {\parallel x\parallel }_0\le k,\hfill \end{array}\right.\end{array}$$

(11)

and

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle \underset{x}{min}{\parallel z-\mathcal{E}Ax\parallel }_2^2}\hfill \\ {\parallel x\parallel }_0\le k.\hfill \end{array}\right.\end{array}$$

(12)

Note that z in (11) and (12) are low-dimensional observations of x obtained from $\mathcal{E}Ax$ and $\mathcal{E}$ in (11) and (12) are the same. To verify that (11) yields performance similar (in terms of probability) to the the sparse recovery obtained using (12), experiments were performed on the sparsifying dictionaries ($D_{wavelet}$ and $D_{KSVD}$) and using the compressive sampling matched pursuit (CoSaMP) algorithm [32] for the recovery of sparse vectors. Figure 3 and Figure 4 illustrate CS recovery performance using plots indicating the probability of successfully recovering a sparse vector versus the CS ratio. The figures compare the performance of our approach (11) versus the benchmark (12) using Gaussian and Bernoulli sensing matrices at various levels of sparsity. The horizontal and vertical axes, respectively, indicate the CS ratio (i.e., ${\displaystyle \frac{m}{n}}\times 100\%$, where m is the number of rows of $\mathcal{E}$ and $n=1024$) and the probability of successfully recovering a sparse vector. We claim that the true sparse vector x can be recovered as long as estimate $\widehat{x}$ satisfies $\parallel \widehat{x}{-x\parallel }_1<1024\times {10}^{-2}$.

5.2. MRI Image Recovery

Significant research effort based on the sparse representation has been directed towards finding ways to accelerate MR imaging. At the heart of these sparse reconstructions assumes MR images can be under-sampled in such a way that data collection time can be dramatically reduced while maintaining image quality. Here, we report experimental results that demonstrate our proposed MRI recovery algorithm on data from the publicly available fastMRI dataset [31]. Our aim is not to develop an algorithm to achieve the state-of-the-art performance for MRI image recovery (baseline deep neural network methods such as variational networks and U-nets [33,34]). The experiment is designed to demonstrate that using the proposed matrix factorization method can be a promising approach for recovering under-sampled MRI images. In this application, one of the factors in our proposed decomposition is adapted to the acceleration mask (restricted by physically prescribed sensing mechanism), which otherwise does not allow for direct RIP-based recovery guarantees in the ${\ell }_1$-synthesis approach to CS. A recent survey of applying the compressed sensing technique to the MRI image recovery [35] presents a variety of techniques to improve the quality of recovered images. In this study, we compared the performance of our method to that based on the TV approach, which promotes sparsity of spacial gradients by approximates a solution to (see [36,37,38]):

$$\begin{array}{c}\hfill \underset{Z\in {\mathbb{R}}^{n_1\times n_2}}{argmin}\parallel Y-\mathcal{R}{F}_{1}Z{F}_2{\parallel }_F^2+\lambda {\displaystyle \sum _{i,j}}{\parallel \nabla Z(i,j)\parallel }_2.\end{array}$$

(13)

Under-sampling masks (restrained by the instrumental conditions) for 4-fold and 8-fold acceleration were applied to full k-space data according to (5). The masks sample 25%, respectively 12.5%, of horizontal scan lines including a fully sampled center region—corresponding to low spatial frequencies—containing 8%, respectively 4%, of scan lines, as shown in Figure 5. We selected CDF 9/7 wavelets with three levels as the sparsifying dictionaries. The lowpass wavelet band and the central areas of the acceleration masks may not be perfectly aligned, leading to potential artifacts in the reconstruction as a result of cross-band interference between the lowpass and highpass wavelet bands.

In the provided images Figure 6, Figure 7, Figure 8 and Figure 9, we can see a visual comparison between the reconstruction outcomes generated by our algorithm and the total-variation (TV) method. These comparisons are based on 4-fold and 8-fold under-sampled Fourier measurements of knee and brain scans from the fastMRI dataset. Parameters were chosen for best performance from a set of values as $\rho =0.02$, $\nu =0.00016$, $\mu =0.00024$, as well as $\lambda =5\times {10}^{-5}$ in (13), and $\gamma =0.0035$ in (10). MRI data: https://fastmri.med.nyu.edu/ (accessed date: 1 July 2022). The parameters were established through empirical testing and adjustment (to obtain the source codes for our algorithm, can contact [email protected] or [email protected]). Images (displayed with 90 degree rotation and cut to quadratic center region) taken from coronal proton density knee dataset (image size 372 × 640, z-slice 20 of file1000031 (knee A) and file1000071 (knee B)), and brain dataset (image size 320 × 640, z-slice 2 of file_brain_AXT1_201_6002786 (brain A) and file_brain_AXT1_201_6002740 (brain B)). The parameters were established through empirical testing and adjustment.

Table 1. Performance comparison for images in Figure 6, Figure 7, Figure 8 and Figure 9, for 4-fold and 8-fold acceleration.

	Proposed				TV Method
	PSNR		SSIM		PSNR		SSIM
	4×	8×	4×	8×	4×	8×	4×	8×
Knee A	28.5	26.4	0.657	0.570	28.2	24.5	0.615	0.538
Knee B	28.4	26.1	0.635	0.535	27.8	24.1	0.625	0.541
Brain A	27.7	24.1	0.642	0.524	26.7	22.3	0.637	0.493
Brain B	27.7	24.6	0.668	0.494	26.8	22.7	0.644	0.474

summarizes a comparison of PSNR and SSIM values for the knee and brain images depicted in the figures. Our average PSNR gains for 4-fold under-sampled data are 0.5 dB on the knee and 0.9 dB on the brain images, while our average gains for the 8-fold under-sampled images are about 2 dB. Our SSIM values are superior to that of the TV-method in all experiments.

6. Conclusions

This article describes a novel approach to CS involving the construction a sensing matrix for a prescribed sparsifying dictionary, such that their composition has the RIP. This property theoretically guarantees the accurate recovery of a fixed number of sparse coefficients. While previous theoretical work has demonstrated RIP assurances only for orthonormal matrices, our study extends these insights to include general over-complete dictionaries. We assert that our sensing matrix can recover sparse coefficients from an over-complete dictionary as effectively as a random sensing matrix for an orthonormal matrix. This advancement is significant, as many practical applications that employ compressive sensing rely on over-complete dictionaries. We further apply the factorization approach to accelerated MR imaging—deriving an embedding for under-sampled k-space data that facilitates RIP recovery guarantees. In this application one of the factors in our proposed decomposition is adapted to the physically prescribed sensing mechanism, which otherwise does not allow for direct RIP based recovery guarantees in the ${\ell }_1$-synthesis approach to CS.

In addition to the achievements mentioned, it is notable to recognize that the analysis could be further improved if the factorization structure in G has a specific arrangement. A structured sensing matrix can significantly enhance storage and computational efficiency in large systems. A reviewer proposes a promising avenue for future exploration for MRI acceleration. Traditional compressive sensing (CS) recovery techniques and modern neural network methods generate different artifacts. While the artifacts from conventional algorithms can usually be traced back to the algorithm and the regularization penalty, those produced by neural networks tend to be more unpredictable. Different neural networks tend to produce distinct model-shifting artifacts. Additionally, the unpredictable nature of the artifacts generated by fixed neural networks presents a significant challenge. Experienced clinicians can quickly identify the artifacts associated with compressive sensing (CS) recovery methods, which helps reduce false positives. However, neural network methods offer a significant advantage, achieving better recovery accuracy and inference speed than traditional techniques. In line with the reviewer’s suggestions, we recommend developing a hybrid approach that combines the strengths and weaknesses of both neural networks and CS techniques to optimize overall performance.