I Introduction
Magnetic resonance (MR) imaging (MRI) provides an indispensable imaging method without ionizing radiation in contemporary clinical applications. However, the scanning speed of MR imaging technique is limited. Both compressed sensing (CS) and parallel imaging (PI) techniques have been applied to reduce the MRI scanning time.
A variety of kspace undersampling patterns has been applied to reduce the amount of collected data, such as onedimensional (1D) uniform undersampling (1DUU), 1D Gaussian random undersampling (1DGU), twodimensional (2D) Poissondisc undersampling (2DPU), 2D Gaussian random undersampling (2DGU), and radial undersampling. According to the CS theory [1, 17], it is feasible to reconstruct the MR images from highly undersampled measurements [17], since the MR images exhibit a sparsity in the wavelet transform domain and spatial finite differences. CSMRI methods [17, 34] have been proposed to solve reconstruction problems with regularization terms of the total variation (TV) and L1 norm of wavelet coefficients. Qu et al. developed the patchbased directional wavelet (PBDW and PBDWS) [24, 21] and graphbased redundant wavelet transform (GBRWT) [13] to improve the reconstruction performance. In addition to fixed sparse transform methods, adaptive sparse representationbased reconstruction methods have been proposed, such as dictionary learningbased MRI (DLMRI) [25] and transform learningbased MRI (TLMRI) [26, 32, 36, 27] algorithms. These algorithms have been proven to attain a high reconstruction performance.
Certain methods have been proposed to exploit the nonlocal selfsimilarity (NSS) of image patches to improve the image quality, such as the nonlocal means (NLmeans) [19], blockmatching 3D denoising (BM3D) [3, 9], patchbased nonlocal operator (PANO) [23], and nonlocal lowrank (NLR)CS [5] methods. The BM3D method exploits the NSS of image patches via their grouping for image denoising purposes [3, 9], which has been applied in MRI reconstruction. Qu et al. [23] exploited the NSS of image patches and established a PANO to reduce the reconstruction error. Dong et al. [5] developed the very promising NLRCS model employing NLR regularization of similar image patches constructed though block matching (BM).
Parallel MRI (PMRI) is also a wellknown technique to accelerate MRI, which is often been combined with the CS theory to improve the reconstruction performance. Sensitivity encoding (SENSE) [22]
explicitly utilizes sensitivity information. However, the performance of these methods is limited due to the difficulty in the accurate measurement of sensitivity information in practical applications. An iterative selfconsistent parallel imaging reconstruction using eigenvector maps (ESPIRiT) model has been proposed
[29]. In contrast to the SENSE model, ESPIRiT model estimates multiple sets of coil sensitivity. We
[6] introduced the Lp pseudonorm joint TV regularization term into the ESPIRiT scheme to improve the reconstruction performance.Another type of PMRI reconstruction method avoids the difficulty in sensitivity estimation by implicitly considering sensitivity information, such as the generalized autocalibrating partially parallel acquisitions (GRAPPA) [11, 2] and iterative selfconsistent parallel imaging reconstruction (SPIRiT) [16, 18, 20, 8, 31, 7] methods, which relies on the kspace local kernel calibration. SPIRiT includes two reconstruction schemes in both the image [18, 20, 31] and kspace [16, 18, 8] domains. An L1SPIRiT scheme was has been obtained by combining the regularization term of the joint L1 norm (JL1) in the wavelet domain with the kspace domainbased SPIRiT model and solved with the projection over convex sets (POCS) algorithm [16, 8]. Duan et al. [7] applied the fast iterative shrinkage thresholding algorithm (FISTA) to solve the SPIRiT PMRI reconstruction problem in the kspace domain with the joint TV (JTV) regularization term. Weller et al. [31] adopted the alternating direction method of multipliers (ADMM) technique to solve the SPIRiT PMRI reconstruction problem in the image domain with the JTV regularization term. Most recently, the STDLRSPIRiT scheme has been established [37] by combining the SPIRiT model with the simultaneous twodirectional lowrankness (STDLR) in kspace to realize improved reconstruction. In fact, STDLRSPIRiT has exploited the local lowrank (LR) prior rather than the LR feature based on nonlocal image structures.
Recent methods have exploited the NSS of images by imposing the group sparsity or the lowrankness of nonlocal similar patches to improve the reconstruction quality. Among them, the NLRbased methods have achieved an excellent performance [5, 33]. In addition, there have been a variety of improved SPIRiTbased methods for PMRI reconstruction. However, to the best of our knowledge, no SPIRiTbased algorithm has applied lowrankness of nonlocal similar patches. In this paper, we propose an NLRSPIRiT scheme incorporating NLR regularization into the SPIRiT model. The NLRSPIRiT model fully utilizes both the NSS in MR images and the calibration consistency in the kspace domain to improve the reconstruction performance. By employing the Nash equilibrium (NE) formulation [4, 35], we reformulate the NLRSPIRiT model into a twoobjective optimization problem including a rank minimization problem and a leastsquares (LS) problem. We adopt the weighted nuclear norm (WNN) instead of the nuclear norm (NN) as a surrogate of the rank, which yields a more efficient method to solve the rank minimization problem. The LS problem is efficiently solved with the ADMM technique. The experimental results demonstrate the superior performance of the NLRSPIRiT model over stateoftheart methods in terms of three objective metrics and visual comparison.
We organize the rest of this article as follows: in Section II, we review the SPIRiT model and SPIRiTbased algorithms. In Section III, we describe the NLRSPIRiT model, which incorporates NLR regularization of similar patches into the SPIRiT model. We adopt the WNN as a surrogate of the rank and employ the NE formulation and ADMM technique to efficiently solve the NLRSPIRiT model. Section IV presents the experimental results and analysis. Finally, we provide the conclusion of this paper in Section V.
Ii Related work
Suppose represents a multicoil image stacked in columns, denotes the undersampled kspace data of the multicoil image stacked in columns, is a undersampled operator selecting only the acquired kspace data from the entire kspace grid,
is the discrete 2D Fourier transform matrix,
is the Fourier operator applied on single coil image, , , represents the undersampled encoding matrix, is a identity matrix, denotes Kronecker product, , and is the total number of coils. The undersampled kspace data of multicoil images are thus given by:(1) 
The SPIRiT calibration consistency equation in the image domain is expressed as:
(2) 
where is the image domainbased SPIRiT operator, acquired from autocalibration signal (ACS) lines (as shown in Fig. 1(a)).
Murphy et al. [20] proposed the image domainbased L1SPIRiT model by introducing the joint sparsity to improve the reconstruction quality. The image domainbased L1SPIRiT minimization problem is as follows:
(3) 
where is the coilbycoil wavelet transform, and denotes joint L1 norm. A simple and efficient POCS algorithm is developed to solve problem (3).
The JL1 and JTV regularization terms are [31] combined with the image domainbased SPIRiT model, and the corresponding optimization problem is given by:
(4)  
where denotes the JTV regularization term. The above optimization problem (4) can be solved with the ADMM technique.
Next, we introduce the kspace domainbased SPIRiT model. The SPIRiT calibration consistency equation in the kspace domain is expressed as:
(5) 
where denotes the kspace data of the multicoil image , and is the kspace domain SPIRiT operator, which convolves a series of calibration kernels acquired from the ACS lines with the entire kspace.
The undersampled kspace data of multicoil images are determined as:
(6) 
Vasanawala et al. proposed the kspacebased L1SPIRT model by introducing the joint sparsity to enhance the reconstruction quality [16]. The L1SPIRiT model in kspace can be expressed as the following minimization problem:
(7) 
Problem (7) is solved with the POCS algorithm [16]. We also discussed the solution to problem (7), and proposed an algorithm with a higher reconstruction quality than that obtained with the POCS algorithm [8].
We introduce the JTV regularization term into the kspacebased SPIRiT model to increase the reconstruction quality [7], and the corresponding optimization problem is given by:
(8)  
By combining the kspacebased SPIRiT model with the STDLR regularization term, Zhang et al. proposed the following STDLRSPIRiT model [37]:
(9)  
where is a fatter matrix comprising stacked Hankel matrices in the horizontal direction, is a fatter matrix comprising stacked Hankel matrices in the vertical direction, is the Hadamard product, is an operator to convert into a block Hankel matrix, and and are the weights corresponding to the Fourier transform coefficients of sparse transform filters in the vertical and horizontal directions, respectively. Problem (9) can be solved with the ADMM technique.
Iii The proposed algorithm
Iiia Problem formulation
Past works have verified that the NSS prior information facilitated sparsity enforcement, resulting in reconstruction quality enhancement. To further improve the PMRI reconstruction quality, we incorporate the NLR regularization term into the SPIRiT model to fully utilize the NSS of PMRI.
Suppose denotes the coil image of , and a single coil image is divided into overlapping patches of size . In each reference patch, we search for similar patches within a local window (e.g., 4040), and choose similar patches closest to the reference patch in terms of the Euclidean distance. The selected similar patches are vectorized and grouped to construct a similar patch group matrix , where is a BM operator. The step of obtaining NSS prior information is depicted in Fig. 1(c). Therefore, the PMRI reconstruction problem based on the NLR regularization term and the SPIRiT model can be formulated as the following problem:
(10)  
where denotes the rank of the matrix , and are tuning parameters, which are applied to balance the data fidelity, the calibration consistency, and the NLR regularization term.
IiiB Problem solution
Problem (10) is a largescale nonconvex optimization problem that is difficult to solve. Authors [4, 35] have proposed employing an NE formulation. With the application of the NE formulation [4, 35], we convert the reconstruction (10) into a twoobjective optimization problem:
(11) 
(12)  
where denotes LR approximation of the patch group matrix , and is determined as follows:
(13) 
where , the adjoint operator of , places back the denoised patches at their original positions (as shown in Fig. 1(c)). is a diagonal matrix with each diagonal element equal to the time of the corresponding pixel belonging to the overlapping patches throughout .
IiiB1 Minimization with respect to
Generally, the rank penalty objective optimization problem of
is a nondeterministic polynomial time (NP)hard problem. Thus, given that the WNN
[12, 15] may yield a better rank approximation than the NN, we adopt the WNN as a convex surrogate of the rank. The WNN of can be written as:(14) 
With the use of the WNN as a surrogate of the rank, let , problem (11) can be rewritten as follows:
(15) 
Problem (15) is a weighted nuclear norm minimization (WNNM) [12] problem. Let
be the full singular value decomposition (SVD) of
, , is the singular value of , and . Hence, the optimal solution to (15) is , , where can be calculated as:(16) 
where is the soft threshold operator, . And the weight can be calculated as:
(17) 
where is a constant, is to avoid dividing by zero, and then the initial can be estimated as:
(18) 
IiiB2 Image reconstruction
After calculating , the whole image can be reconstructed by solving problem (12). By introducing an auxiliary variable and corresponding Lagrange multiplier , the LS problem (12) can be converted to the following subproblems via the ADMM technique:
(19) 
(20)  
(21) 
The solution to subproblem (19) with respect to is given by:
(22) 
Let , and can thus be obtained via the direct inversion of each block in matrix , as shown in Fig. 1(b). Then, we have the update of as . is a diagonal matrix, and subproblem (20) with respect to yields the following closedform solution:
(23) 
Since subproblems (11), (19), and (20) are efficiently solved, the Ccoil image is obtained, and then combined into a single magnitude image by using the square root of sum of squares (SOS):
(24) 
Finally, we obtain the NLR regularizationbased SPIRiT (NLRSPIRiT) PMRI reconstruction algorithm, as expressed in Algorithm 1. And the schematic illustration of NLRSPIRiT is described in Fig. 1. Algorithm 1 is terminated when the relative error (RE) falls below the tolerance . In Algorithm 1, the BM step is performed every iterations to reduce the computational complexity, as shown in step 4.
The optimization problem (10) can also be solved by using the variable splitting and ADMM techniques. The corresponding algorithm is called ADMMNLRSPIRiT. Please refer to the appendix for details.
Iv Experimental results
Iva Experimental setup
In our experiments, we compared the proposed NLRSPIRiT model to stateoftheart algorithms to solve the PMRI reconstruction problems such as: JTVSPIRiT [31], LPJTVESPIRiT [6], STDLRSPIRiT [37], and NLRSPIRiTbaseline (a variant of the NLRSPIRiT model with the standard NN). All the considered algorithms were implemented in MATLAB.
To validate the performances of all the considered algorithms, we conduct experiments on four fully sampled in vivo human datasets. The first fully sampled brain dataset [28] (dataset 1, as shown in Fig. 3LABEL:sub@fig3a) was acquired on a clinical MR scanner Discovery MR750 using a 12channel headneck spine coil and a 3D T1weighted gradient echo sequence (matrix size = , TR/TE = 7.4/3.1 ms, field of view (FOV) = mm, slice thickness = 1 mm). We chose a single slice of each multislice dataset in our experiment. The second fully sampled knee dataset [10] (dataset 2, as shown in Fig. 4LABEL:sub@fig4a) was acquired on a GE scanner using an 8channel HD knee coil and a 3D FSE CUBE sequence (matrix size = , TR/TE = 1550/25 ms, FOV = mm, slice thickness = 0.6 mm). We chose a single slice of each multislice dataset in our experiment. The third fully sampled brain dataset [14] (dataset 3, as shown in Fig. 5LABEL:sub@fig5a) was acquired on a 3T SIEMENS Trio system and an 8channel head array coil (matrix size = , TR/TE = 2530/3.45 ms, slice thickness = 1.33 mm, FOV = mm). The fourth fully sampled brain dataset [37] (dataset 4, as shown in Fig. 6LABEL:sub@fig6a) was acquired on a 3T SIEMENS Trio wholebody scanner using a 32coil head coil and the 2D T2weighted turbo spin echo sequence (matrix size = , TR/TE = 6100/99 ms, FOV = mm, slice thickness = 3 mm). The 32channel kspace data was compressed to four virtual coils.
To validate the considered algorithms, fully sampled datasets were subjected to retrospective undersampling with the above different undersampling patterns and acceleration factors (AFs) (including ACS lines). In our experiments, we employed two undersampling patterns, 2DPU and 1DGU. It should be noted that the 2DPU patterns always include 2424 ACS lines, and the 1DGU pattern includes 20 ACS lines. All experiments were simulated on a workstation machine equipped with an Intel Core i910900X @ 3.7 GHz processor, 64 G of RAM memory and a Windows 10 operating system (64 bit).
Three commonly adopted objective metrics were employed to evaluate the quality of the reconstructed images, namely, the signaltonoise ratio (SNR)
[7], highfrequency error norm (HFEN) [25], and structural similarity index measure (SSIM) [30]. It should be noted that the calculation of all metrics is limited to the region of interest (ROI). High SNR and SSIM values or low HFEN values indicate more accurate reconstruction. In regard to reference image and the reconstructed image , the SNR is defined as:(25) 
where denotes the mean square error between and , and
is the variance in
.The HFEN is expressed as:
(26) 
where is a Laplacian Gaussian filter operator used to determine the image edges.
The SSIM is calculated as:
(27) 
where and are the means of and , respectively, and are the variances of and , respectively, represents the covariance of and , and are constant.
IvB Parameter settings
A kernel size of was adopted in the SPIRiTbased algorithms in the following experiments. The parameters of JTVSPIRiT, LPJTVESPIRiT, and STDLRSPIRiT were manually tuned for SNR optimal. In regard to NLRSPIRiTbaseline and NLRSPIRiT, image patches were used, and the total number of similar patches was 43 . A reference image patch was extracted every 5 pixels along the horizontal and vertical directions to reduce the computational complexity. The main parameters were defined as follows: , , and . These settings remained fixed in all experiments.
In addition, to obtain a better parallel image reconstruction performance, the parameter and were adjusted to optimize the SNR and HFEN values. We run the proposed algorithm on the same undersampled data considering the ranges of and . We chose and based on the highest SNR and smallest HFEN values. For example, at and , these two optimal indicators (SNR and HFEN) of dataset 1 were roughly optimized.
Algorithms  Metrics  Acceleration Factor (AF)  

3  4  5  6  7  
JTVSPIRiT [31]  SNR  24.35  22.85  21.55  20.38  19.56 
HFEN  0.0938  0.1197  0.1505  0.1816  0.2048  
SSIM  0.9679  0.9585  0.9483  0.9373  0.9277  
LPJTVESPIRiT [6]  SNR  24.24  22.73  21.50  20.33  19.53 
HFEN  0.0942  0.1203  0.1479  0.1802  0.2008  
SSIM  0.9696  0.9593  0.9341  0.9363  0.9264  
STDLRSPIRiT [37]  SNR  24.54  22.91  21.45  20.22  19.34 
HFEN  0.0898  0.1141  0.1438  0.1715  0.1943  
SSIM  0.9689  0.9577  0.9452  0.9328  0.9220  
NLRSPIRiTbaseline  SNR  24.79  23.37  22.16  21.16  20.36 
HFEN  0.0898  0.1122  0.1381  0.1605  0.1804  
SSIM  0.9712  0.9626  0.9532  0.9436  0.9344  
NLRSPIRiT  SNR  25.37  24.16  22.91  21.90  21.05 
HFEN  0.0810  0.0995  0.1244  0.1509  0.1744  
SSIM  0.9749  0.9690  0.9612  0.9516  0.9429 
Algorithms  Metrics  Acceleration Factor (AF)  

3  4  5  6  7  
JTVSPIRiT [31]  SNR  18.75  17.53  16.85  16.23  15.76 
HFEN  0.2507  0.3022  0.3272  0.3568  0.3855  
SSIM  0.9161  0.8953  0.8831  0.8739  0.8638  
LPJTVESPIRiT [6]  SNR  19.20  17.97  17.41  16.75  16.20 
HFEN  0.2110  0.2552  0.2801  0.3062  0.3349  
SSIM  0.9233  0.9014  0.8307  0.8781  0.8663  
STDLRSPIRiT [37]  SNR  19.34  17.83  17.02  16.07  15.24 
HFEN  0.2245  0.2810  0.3084  0.3483  0.3886  
SSIM  0.9262  0.9026  0.8878  0.8714  0.8524  
NLRSPIRiTbaseline  SNR  19.13  18.42  17.75  17.20  16.67 
HFEN  0.2146  0.2427  0.2655  0.2909  0.3143  
SSIM  0.9201  0.9094  0.8984  0.8889  0.8781  
NLRSPIRiT  SNR  20.24  19.08  18.33  17.75  17.19 
HFEN  0.1860  0.2184  0.2421  0.2667  0.2921  
SSIM  0.9352  0.9172  0.9049  0.8947  0.8852 
Algorithms  Metrics  Acceleration Factor (AF)  

3  4  5  6  7  
JTVSPIRiT [31]  SNR  27.36  24.73  22.65  20.55  19.22 
HFEN  0.1362  0.1836  0.2322  0.2994  0.3382  
SSIM  0.9756  0.9619  0.9463  0.9273  0.9112  
LPJTVESPIRiT [6]  SNR  27.46  25.06  23.11  21.22  20.02 
HFEN  0.1322  0.1734  0.2162  0.2743  0.3096  
SSIM  0.9810  0.9704  0.9585  0.943  0.9291  
STDLRSPIRiT [37]  SNR  27.53  25.11  22.99  20.94  19.48 
HFEN  0.1333  0.1765  0.2245  0.2882  0.3316  
SSIM  0.9807  0.9702  0.9576  0.9409  0.9250  
NLRSPIRiTbaseline  SNR  27.74  25.36  23.80  21.80  20.45 
HFEN  0.1294  0.1685  0.2022  0.2586  0.2948  
SSIM  0.9819  0.9727  0.9644  0.9497  0.9381  
NLRSPIRiT  SNR  28.11  26.11  24.31  22.27  20.73 
HFEN  0.1223  0.1544  0.1909  0.2446  0.2849  
SSIM  0.9836  0.9761  0.9673  0.9529  0.9406 
Algorithms  SNR  HFEN  SSIM 

JTVSPIRiT [31]  23.26  0.0688  0.9901 
LPJTVESPIRiT [6]  23.77  0.0578  0.9906 
STDLRSPIRiT [37]  24.46  0.0665  0.9919 
NLRSPIRiTbaseline  24.58  0.0608  0.9925 
NLRSPIRiT  25.51  0.0536  0.9935 
IvC Convergence analysis
To reflect the convergence of the proposed NLRSPIRiT algorithm, Fig. 2 shows SNR, HFEN, and RE versus the iteration number during the reconstruction of dataset 1 based on the 2DPU pattern . We observe that NLRSPIRiT approximately reaches the maximum SNR and the minimum HFEN when RE falls below 1e4. At this time, NLRSPIRiT converges and then is terminated when RE falls below 1e4.
IvD Comparison to previous works
In the following experiments, JTVSPIRiT [31], LPJTVESPIRiT [6], STDLRSPIRiT [37], NLRSPIRiTbaseline, and NLRSPIRiT were compared.
IvD1 Comparison results of 2DPU patterns
To facilitate the assessment of the image quality (commonly a subjective parameter), Figs. 35 show the images reconstructed via all considered algorithms based on the 2DPU patterns , and their corresponding error maps. Fig. 3 shows that the reconstructed images via JTVSPIRiT obvious artifacts. LPJTVESPIRiT, STDLRSPIRiT and NLRSPIRiTbaseline mitigate these artifacts to a certain extent. The proposed NLRSPIRiT algorithm effectively removes image artifacts and preserves the details. Hence, NLRSPIRiT achieves the best visual quality among all competing algorithms. Fig. 4 shows that the reconstructed images via JTVSPIRiT and NLRSPIRiTbaseline exhibit apparent artifacts. LPJTVESPIRiT and STDLRSPIRiT mitigate artifacts to a certain extent. The proposed NLRSPIRiT algorithm further removes these artifacts and succeeds in reconstructing certain details accurately, such as knee joint, whereas artifact areas remain in images more obviously reconstructed via other considered algorithms. NLRSPIRiT produces the visually best reconstructed image based on visual comparison. Fig. 5 shows that the reconstructed images via JTVSPIRiT and STDLRSPIRiT exhibit apparent artifacts. NLRSPIRiTbaseline and STDLRSPIRiT mitigate these artifacts effectively. NLRSPIRiT obtain fewest errors and produces the visually best reconstructed image based on visual comparison. In summary, the proposed NLRSPIRiT model achieves the best visual performance.
To quantitatively evaluate the reconstruction performance of all considered algorithms based on 2DPU patterns with different AF values, comparison results of the SNR, HFEN, and SSIM values of the reconstructed images via the considered algorithms are summarized in Tables IIII. The best metric in each case is marked in bold.
As indicated in Tables IIII, one can see that the proposed NLRSPIRiT method produces the best metrics for all datasets and AFs. For dataset 1, NLRSPIRiT achieves SNR improvements of 1.34, 1.41, 1.39, and 0.71 dB on average over JTVSPIRiT, LPJTVESPIRiT, STDLRSPIRiT, and NLRSPIRiTbaseline, respectively. In regard to dataset 2, NLRSPIRiT achieves SNR improvements of 1.50, 1.01, 1.42, and 0.68 dB, respectively, on average. Regarding dataset 3, NLRSPIRiT achieves SNR improvements of 1.40, 0.93, 1.09, and 0.48 dB, respectively, on average.
NLRSPIRiTbaseline attains a better reconstruction performance than that of JTVSPIRiT, LPJTVESPIRiT, and STDLRSPIRiT except for a case (dataset 2 with ). In other words, the introduction of the NLR regularization term contributes to reconstruction performance improvement of these SPIRiTbased algorithms. Furthermore, NLRSPIRiT improve the reconstruction performance than NLRSPIRiTbaseline. NLRSPIRiT achieves a significant improvement over JTVSPIRiT, LPJTVESPIRiT, and STDLRSPIRiT in terms of the three metrics.
IvD2 Comparison results of 1DGU patterns
We choose the same data (dataset 4) and sampling pattern (1DGU pattern at a sampling rate of 0.34, which derived from Fig. 7(f) in 22) for comparison with the STDLRSPIRiT [37] algorithm. The comparison of the PMRI reconstruction performance for dataset 4 under the 1DGU pattern is shown in Fig. 6LABEL:sub@fig6f.
As shown in Fig. 6, all the considered algorithms can produce a good reconstruction. And NLRSPIRiT achieves the least artifacts inside the image and the least errors at the edge. It is shown that the proposed NLRSPIRiT algorithm yields the highest visual quality among all considered algorithms.
IvD3 Statistical analysis
We also provide violin plots for the relative metric differences between NLRSPIRiT and the compared algorithms in Table IIV. As shown in Fig. 7, NLRSPIRiT significantly outperforms JTVSPIRiT, LPJTVESPIRiT, and STDLRSPIRiT in terms of SNR, HFEN, and SSIM. In addition, we can also see that NLRSPIRiT has slightly better performance than NLRSPIRiTbaseline.
IvD4 Comparison of runtime and memory demand
Next, we compare the runtime and memory demands of JTVSPIRiT, LPJTVESPIRiT, STDLRSPIRiT, and NLRSPIRiT. Table V summarizes the runtime and memory demands of the considered algorithms in PMRI reconstruction for dataset 1 based on the 2DPU pattern (). As indicated in Table V, JTVSPIRiT has advantages in runtime and memory demands. Compared to JTVSPIRiT, LPJTVESPIRiT consumes fewer memory but requires more runtime. STDLRSPIRiT is the most advanced SPIRiTbased method at hand, which has the highest computational complexity, requires the longest runtime, and the maximum memory. Furthermore, NLRSPIRiT requires 0.49 GB of RAM memory for reconstruction, only approximately one92th of the memory demand of STDLRSPIRiT. In other words, NLRSPIRiT outperforms STDLRSPIRiT in terms of memory demand and runtime.
The runtime of NLRSPIRiT is , where is the outer iteration number, is the runtime of the blockmatching operation performed every iterations, is the runtime of one LR approximation, is the runtime to update , and once, and is the runtime of algorithm initialization and preprocessing. The very high computational complexity of the LR approximation and blockmatching steps result in large and values. Therefore, we will optimize the LR approximation step and blockmatching steps with a graphics processing unit (GPU) to efficiently accelerate the proposed NLRSPIRiT algorithm in the future.
IvE Proposed NE formulation versus ADMM formulation
According to the derivations of the NLRSPIRiT model and ADMMNLRSPIRiT algorithm (see the appendix), we observe that these two algorithms attain almost the same calculation time for oneiteration.
We compared the reconstruction results obtained with NLRSPIRiT and ADMMNLRSPIRiT for dataset 1 based on the 2DPU pattern . As shown in Fig. 8, NLRSPIRiT attains faster convergence and obtains a slightly higher SNR value than ADMMNLRSPIRiT. And we determine that they require almost the same reconstruction time every iteration. In addition, NLRSPIRiT contains two fewer parameters than does ADMMNLRSPIRiT, and therefore it is easier to tune the parameters of NLRSPIRiT. In summary, it is concluded that the proposed NLRSPIRiT algorithm exhibits obvious advantages in practical applications.
IvF Further discussion about parameter settings
We found that there exist similar parameters for reconstructing the MR images of same kind, so we also offer a parameter selection method for practical application. First, we train the parameters from 15 retrospective undersampling data sets with the 2DPU pattern (AF = 5). We ran NLRSPIRiT on the training undersampling data sets on the ranges of and , and obtained SNRs for different parameter settings. The average SNR relative difference between the maximum SNR and the SNR for each parameter setting is calculated. As shown in Fig. 9(a), the average SNR relative differences are minimized at and . Second, we ran NLRSPIRiT on the 32 retrospective validation undersampling data sets on the ranges of and , and obtained the average SNR relative differences for each parameter setting. As shown in Fig. 9(b), the average SNR relative differences are also minimized around and . Therefore, for this kind of data sets, the optimal parameter setting was determined to be and .
V Conclusions
In this paper, we propose the NLRSPIRiT model, which incorporates the NLR regularization into the SPIRiT model. The NLRSPIRiT model fully utilizes both the NSS in MR images and the calibration consistency in the kspace domain. We adopt the WNN instead of the NN as a surrogate of the rank, and employ the NE formulation and the ADMM technique to efficiently solve the NLRSPIRiT model. Experimental results considering four vivo datasets indicate that the proposed NLRSPIRiT algorithm almost achieves a superior performance in terms of three objective metrics and visual perception over stateoftheart methods. In addition, we propose a parameter setting method for practical application, which selects a set of nearoptimal parameters for the same kind of MR images. We will optimize the most timeconsuming BM and LR approximation steps with a GPU to efficiently accelerate the NLRSPIRiT algorithm in the future. The proposed NLRSPIRiT algorithm is very promising in regard to PMRI applications.
Acknowledgment
The authors would like to acknowledge Michael Lustig, Daniel S. Weller, Xinlin Zhang, Kevin Epperson, FaHsuan Lin, Robreto Souza, Saiprasad Ravishankar, Bihan Wen, and Weisheng Dong, for making their code or in vivo data sets publicly available. The authors would like to thank the anonymous reviewers for their comments and suggestions.
Appendix A Alternative algorithm
The PMRI reconstruction problem based on the NLR regularization and the SPIRiT model can be rewritten as follows:
(28)  
Problem (28) can be solved with the variable splitting (VS) and ADMM techniques. First of all, by introducing auxiliary variables , , and , in addition to corresponding Lagrange multipliers , , , and , respectively, problem (28) is decomposed into the following subproblems via the ADMM method:
(29) 
(30) 
(31)  
(32)  
(33) 
(34) 
(35) 
As mentioned in Section III of the manuscript, the subproblems (29) and (30) are efficiently solved with respect to and , respectively. Subproblem (31) with respect to yields the following closedform solution:
(36) 
Subproblem (32) with respect to generates a closedform solution, which is expressed as follows:
(37) 
References
 [1] (2006) Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Transactions on Information Theory 52 (2), pp. 489–509. External Links: Document Cited by: §I.
 [2] (2012) Nonlinear GRAPPA: a kernel approach to parallel MRI reconstruction. Magnetic Resonance in Medicine 68 (3), pp. 730–740. External Links: Document Cited by: §I.
 [3] (2007) Image denoising by sparse 3D transformdomain collaborative filtering. IEEE Transactions on Image Processing 16 (8), pp. 2080–2095. External Links: Document Cited by: §I.
 [4] (201204) BM3D frames and variational image deblurring. IEEE Trans Image Process 21 (4), pp. 1715–1728. External Links: Document Cited by: §I, §IIIB.
 [5] (2014) Compressive sensing via nonlocal lowrank regularization. IEEE Transactions on Image Processing 23 (8), pp. 3618–3632. External Links: Document Cited by: §I, §I.
 [6] (2019) Eigenvectorbased SPIRiT parallel MR imaging reconstruction based on Lp pseudonorm joint total variation. Magnetic Resonance Imaging 58, pp. 108–115. External Links: Document Cited by: §I, Fig. 3, Fig. 4, Fig. 5, Fig. 6, §IVA, §IVD, TABLE I, TABLE II, TABLE III, TABLE IV, TABLE V.
 [7] (2018) Efficient operator splitting algorithm for joint sparsityregularized SPIRiTbased parallel MR imaging reconstruction. Magnetic Resonance Imaging 46, pp. 81–89. External Links: Document Cited by: §I, §II, §II, §IVA.
 [8] (2014) Efficient reconstruction algorithm for parallel magnetic resonance imaging based on selfconsistency. Journal of Tianjin University 47 (5), pp. 414–419. External Links: Document Cited by: §I, §II.
 [9] (2016) Decoupled algorithm for MRI reconstruction using nonlocal block matching model: BM3DMRI. Journal of Mathematical Imaging and Vision 56 (3), pp. 430–440. External Links: Document Cited by: §I.
 [10] (2013) Creation of fully sampled MR data repository for compressed sensing of the knee. Cited by: §IVA.
 [11] (2002) Generalized autocalibrating partially parallel acquisitions (GRAPPA). Magnetic Resonance in Medicine 47 (6), pp. 1202–1210. External Links: Document Cited by: §I.

[12]
(2014)
Weighted nuclear norm minimization with application to image denoising.
In
IEEE Conference on Computer Vision and Pattern Recognition (CVPR)
, Columbus, OH, USA, pp. 2862–2869. External Links: Document Cited by: §IIIB1, §IIIB1.  [13] (2016) Image reconstruction of compressed sensing MRI using graphbased redundant wavelet transform.. Medical Image Analysis 27, pp. 93–104. External Links: Document Cited by: §I.
 [14] Mprage_8ch_slice1.mat. Note: http://www.nmr.mgh.harvard.edu/~fhlin/codes/mprage_8ch_slice1.matAccessed April 4, 2013 Cited by: §IVA.
 [15] (2016) Nonconvex nonsmooth low rank minimization via iteratively reweighted nuclear norm. IEEE Transactions on Image Processing 25 (2), pp. 829–839. External Links: Document Cited by: §IIIB1.
 [16] (2009) L1 SPIRiT: autocalibrating parallel imaging compressed sensing. In International Society for Magnetic Resonance in Medicine (ISMRM), Honolulu, Hawaii, pp. 334. Cited by: §I, §II, §II.
 [17] (2007) Sparse MRI: the application of compressed sensing for rapid MR imaging. Magnetic Resonance in Medicine 58 (6), pp. 1182–1195. External Links: Document Cited by: §I.
 [18] (2010) SPIRiT: iterative selfconsistent parallel imaging reconstruction from arbitrary kspace. Magnetic Resonance in Medicine 64 (2), pp. 457–471. External Links: Document Cited by: §I.
 [19] (2010) Adaptive nonlocal means denoising of MR images with spatially varying noise levels. Journal of Magnetic Resonance Imaging 31 (1), pp. 192–203. External Links: Document Cited by: §I.
 [20] (2012) Fast L1SPIRiT compressed sensing parallel imaging MRI: scalable parallel implementation and clinically feasible runtime. IEEE Transactions on Medical Imaging 31 (6), pp. 1250–1262. External Links: Document Cited by: §I, §II.
 [21] (2013) Magnetic resonance image reconstruction using trained geometric directions in 2D redundant wavelets domain and nonconvex optimization.. Magnetic Resonance Imaging 31 (9), pp. 1611–1622. External Links: Document Cited by: §I.
 [22] (199911) SENSE: sensitivity encoding for fast MRI. Magnetic Resonance in Medicine 42 (5), pp. 952–962. External Links: Document Cited by: §I.
 [23] (2014) Magnetic resonance image reconstruction from undersampled measurements using a patchbased nonlocal operator. Medical Image Analysis 18 (6), pp. 843–856. External Links: Document Cited by: §I.
 [24] (2012) Undersampled MRI reconstruction with patchbased directional wavelets.. Magnetic Resonance Imaging 30 (7), pp. 964–977. External Links: Document Cited by: §I.
 [25] (201105) MR image reconstruction from highly undersampled kspace data by dictionary learning. IEEE Transactions on Medical Imaging 30 (5), pp. 1028–1041. External Links: Document Cited by: §I, §IVA.
 [26] (2015) Efficient blind compressed sensing using sparsifying transforms with convergence guarantees and application to magnetic resonance imaging. SIAM Journal on Imaging Sciences 8 (4), pp. 2519–2557. External Links: Document Cited by: §I.
 [27] (201712) Efficient sum of outer products dictionary learning (SOUPDIL) and its application to inverse problems. IEEE Transactions on Computational Imaging 3 (4), pp. 694–709. External Links: Document Cited by: §I.
 [28] (202009) Dualdomain cascade of Unets for multichannel magnetic resonance image reconstruction. Magnetic Resonance Imaging 71, pp. 140–153. External Links: Document Cited by: §IVA.

[29]
(2014)
ESPIRiTan eigenvalue approach to autocalibrating parallel MRI: where SENSE meets GRAPPA
. Magnetic Resonance in Medicine 71 (3), pp. 990–1001. External Links: Document Cited by: §I.  [30] (200404) Image quality assessment: from error visibility to structural similarity. IEEE Transactions on Image Processing 13 (4), pp. 600–612. External Links: Document Cited by: §IVA.
 [31] (2014) Augmented Lagrangian with variable splitting for faster nonCartesian L1SPIRiT MR image reconstruction. IEEE Transactions on Medical Imaging 33 (2), pp. 351–361. External Links: Document Cited by: §I, §II, Fig. 3, Fig. 4, Fig. 5, Fig. 6, §IVA, §IVD, TABLE I, TABLE II, TABLE III, TABLE IV, TABLE V.
 [32] (2015) FRIST  flipping and rotation invariant sparsifying transform learning and applications. Inverse Problems 33 (7), pp. 074007. External Links: Document Cited by: §I.
 [33] (2020) Image recovery via transform learning and lowrank modeling: the power of complementary regularizers. IEEE Transactions on Image Processing 29, pp. 5310–5323. External Links: Document Cited by: §I.
 [34] (2010) A fast alternating direction method for TVL1L2 signal reconstruction from partial Fourier data. IEEE Journal of Selected Topics in Signal Processing 4 (2), pp. 288–297. External Links: Document Cited by: §I.
 [35] (2014) Motion adaptive patchbased lowrank approach for compressed sensing cardiac cine MRI. IEEE Transactions on Medical Imaging 33 (11), pp. 2069–2085. External Links: Document Cited by: §I, §IIIB.
 [36] (2016) Fast multiclass dictionaries learning with geometrical directions in MRI reconstruction. IEEE Transactions on Biomedical Engineering 63 (9), pp. 1850–1861. External Links: Document Cited by: §I.
 [37] (2020) Image reconstruction with lowrankness and selfconsistency of kspace data in parallel MRI. Medical Image Analysis 63, pp. 101687. External Links: Document Cited by: §I, §II, Fig. 3, Fig. 4, Fig. 5, Fig. 6, §IVA, §IVA, §IVD2, §IVD, TABLE I, TABLE II, TABLE III, TABLE IV, TABLE V.