{'fraction_non_alphanumeric': 0.041700885117255226, 'fraction_numerical': 0.015664446269428475, 'mean_word_length': 4.651082846338948, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accelerated MRI Reconstruction with Separable and Enhanced Low-Rank Hankel Regularization

The combination of the sparse sampling and the low-rank structured matrix reconstruction has shown promising performance, enabling a significant reduction of the magnetic resonance imaging data acquisition time. However, the low-rank structured approaches demand considerable memory consumption and are time-consuming due to a noticeable number of matrix operations performed on the huge-size block Hankel-like matrix. In this work, we proposed a novel framework to utilize the low-rank property but meanwhile to achieve faster reconstructions and promising results. The framework allows us to enforce the low-rankness of Hankel matrices constructing from 1D vectors instead of 2D matrices from 1D vectors and thus avoid the construction of huge block Hankel matrix for 2D k-space matrices. Moreover, under this framework, we can easily incorporate other information, such as the smooth phase of the image and the low-rankness in the parameter dimension, to further improve the image quality. We built and validated two models for parallel and parameter magnetic resonance imaging experiments, respectively. Our retrospective in-vivo results indicate that the proposed approaches enable faster reconstructions than the state-of-the-art approaches, e.g., about 8x faster than STDLRSPIRiT, and faithful removal of undersampling artifacts.

artifacts into the images and may not be practical in some application scenarios [1].
To alleviate the prolonged acquisition time, researchers have made great efforts. Shortening the scan time by acquire only a small subset of the k-space data has emerged as an effective and widely-used way [2], [3]. The sparse sampling technique enables significant reduction of time, but results in undersampling artifacts in the images due to sub-Nyquist sampling. Therefore, reconstruction models with prior information are exploited to enable promising results. Sparsity, lowrank, and sparsity plus low-rank are commonly-used priors. For instance, the sparse MRI methods enforce the sparsity of images in the transform domain, e.g., total variation [4], wavelets [2], contourlets [5], and adaptive sparse transform [6]- [11]. The low-rank prior stems from the linear correlations among multiple MRI images [3], such as, dynamic MRI [12], [13], parameter imaging [14], [15], and high-dimensional MRI [16], [17] and has been utilized.
The low-rank structured matrix prior has been shown to be able to permit promising results [18]- [22]. For instance, simultaneous autocalibrating and k-space estimation [19] utilizes the strong correlation among receiver coils to enable the low-rankness of a block Hankel matrix constructed from the multi-coil k-space data. The low-rank matrix modeling of local k-space neighborhoods (LORAKS) [18] assumes the image of interest has limited support and/or a slowly varying phase, which would lead to a block-Hankel-like matrix being low-rank. The annihilating filter-based low-rank Hankel matrix approach (ALOHA) [20], and Ongie's method [21] consider the sparsity of the signal in the transform domain. Similarly, the simultaneous two-directional low-rankness with SPIRiT (STDLR-SPIRiT) method [22] exploits the low-rank prior and enforces the self-consistency of k-space data providing stateof-the-art performance.
The structured low-rank methods lift the signal into a higher dimensional space to achieve low-rankness. However, the lifting brings some limitations, such as, large memory consumption and lengthy computational time. As for four-coil 256 × 256 data, the dimension of its block Hankel matrix reaches 2116×54756 when the pencil parameters equal 23. Researcher have established methods aiming to tackle the problem. The generalized iterative reweighted annihilating filter algorithm [23] utilized the convolutional structure of the lifted matrix so that it can performs the computations in the original signal instead of the lifted matrix. This method dramatically reduces the complexity and shortens the computational time. Auto-calibrated-LORAKS (AC-LORAKS) [24] exploited the auto-calibration signal (ACS) to reduce the computational time. Though the two methods gain considerable acceleration in time, they take advantage of specific priors that may not available in other low-rank approaches. Thus, fast low-rank Hankel MRI reconstruction is still highly desirable.
In this work, we attempted to establish a framework to permit fast and reliable low-rank reconstructions. A separable low-rank reconstruction method was proposed by enforcing the low-rankness of multiple small Hankel matrices from each row and each column of MRI data, avoiding constructing the high-dimensional low-rank matrix thereby enables much less memory and allows faster computation. In addition, the self-consistence and conjugate symmetric of k-space data are considered to further reduce the reconstruction error. Besides, the parameter-dimensional information will be introduced into formulation so that the original model can not only be applied to non-parameter parallel imaging but also extended to parameter imaging. The proposed method is expected to effectively utilize the low-rank property to yield fast reconstruction and maintain low reconstruction error.
The rest of this paper is organized as follows. We introduce some notations in Section II, and related work in Section III. Section IV presents the proposed model and numerical algorithm for non-parameter and parameter parallel imaging reconstructions, and Section III demonstrates the reconstruction performance. Section VI discusses some factors that may affect the reconstruction error and speed. The conclusions are finally drawn in Section VII.

II. NOTATIONS
For easier reading, we first list some notations used throughout this paper in Table I. III. RELATED WORK The STDLR [22] simultaneously enforced the low-rank constraint of k-space data along both horizontal and vertical directions as follows:

Symbol
Quantity Note: the FE is the abbreviation of frequency encoding, the PE is the abbreviation of phase encoding, the capital letter P denotes the parameter dimension.
where K ∈ C M ×N ×J denotes the targeted k-space data, and Y ∈ C M ×N ×J the acquired multi-coil k-space data with zerofilling in unacquired position. The operator B transfers multicoil data into a cascaded block Hankel matrix; W 2D = and W 2D ⊥ are weighting operators that perform weighting on kspace data of each coil image, whose weights are the Fourier transform of filters in the horizontal and vertical directions; U represents an operator that undersamples data, and vec (·) means arranging a matrix or a tensor to a vector.
The STDLR explored the simultaneous horizontal and vertical directional low-rankness in the k-space data, reducing the image reconstruction errors. Compared to the state-of-the-art method ALOHA, STDLR achieves lower reconstruction errors [22]. However, the size of the block Hankel matrix appears dramatically huge, leading to a long computational time and large memory requirements. Take four-coil parallel imaging data for an example, the STDLR requires 1014.2 s to finish reconstruction while the separable Hankel low-rank method takes only 40.8 s (Fig. 2).

A. Basic Model -SHLR
In this work, we first proposed a separable Hankel low-rank reconstruction method (SHLR) to enforce the low-rankness of Hankel matrices constructing from each row and each column as follows: where the X ∈ C M ×N ×J denotes the desired multi-coil image, F 2D is the 2D Fourier transform for each coil data. The P m and Q n denote the operators that extract m-th row (m-th vector on the first dimension) and n-th column (n-th vector on the second dimension) from each coil data for m = 1, · · · , M and n = 1, · · · , N . Here, we define x row m,j as the m-th row and j-th coil of X, and x col n,j as the n-th column and j-th coil of X. Then, we have P m X = x row m,1 , · · · , x row m,j , · · · , x row m,J ∈ C N ×J and Q n X = x col n,1 , · · · , x col n,j , · · · , x col n,J ∈ C M ×J . Fig. 3 shows a graphical illustration of the notations of constructing Hankel matrices. The operator F 1D denotes the 1D Fourier transform on a vector, W performs weighting on a vector with the weights obtained from applying Fourier transform to 1D sparse transform filter, and H converts a vector into a Hankel matrix. The tilde above the operator means that corresponding operation is performed on each column vector of the matrix, that isHWF 1D P m X = HWF 1D x row m,1 , · · · , HWF 1D x row m,J andHWF 1D Q n X = HWF 1D x col n,1 , · · · , HWF 1D x col n,J .

B. Enhanced Model
As mentioned above, building the block Hankel matrix requires huge memory consumption and thereby increases the computational complexity and time. To mitigate the problem, we proposed to establish multiple small-size Hankel matrices by exploiting the low-rankness of a row or a column vector of the signal. This will significantly reduce memory consumption and make computational speed much faster.
However, constraining the low-rankness of each row and each column separately weakens the utilization between rows and columns, leading to an increase in reconstruction errors, e.g. relative 2 norm error (RLNE). (Fig. 2 (e-f)). In order to improve the reconstruction while maintaining a faster reconstruction speed, we consider introducing other information to strengthen the constraints of the model. The following three kinds of information are exploited.

1) Strengthen the correlation between rows and columns:
The proposed SHLR method makes little use of the correlation between rows and columns. A natural idea for improvement is to strengthen the exploitation of k-space neighborhood information. The SPIRiT is a good choice. The SPIRiT [25] primarily bases on the assumption that each k-space data point is the convolution of multi-coil data of its neighboring kspace, which implies the SPIRiT constraint has the ability to utilize the linear correlation of adjacent rows and columns to reconstruct MRI image. Therefore, the SHLR with SPIRiT constraint (SHLR-S) can be modeled as: where G denotes the SPIRiT operator in image domain.
2) Strengthen the low-rankness within each row (or column): Another way to improve the reconstruction is to strengthen the low-rankness within each row and each column. The smooth phase of MRI images is often used as prior information to be involved in MRI reconstruction [18], [26], [27]. By generating the virtual coils containing the conjugate symmetric, the smooth MRI image phase information can be incorporated into the reconstruction problem [28]. The conjugate symmetry of k-space has been applied into the Hankel low-rank reconstruction [29]. This property can be easily introduced into SHLR reconstruction model as: HvcWF 1D QnX = HWF 1D x col n,1 , · · · , HWF 1D x col n,J , where the superscript † represents the operation of taking the conjugate and flipping the vector along the center. Therefore, we can formulate the SHLR with virtual coil (SHLR-V) as: We introduced both the SPIRiT and virtual coil information into the basis SHLR model aiming to further reduce the reconstruction errors. We can achieve better reconstruction results with the optimization as: As shown in Fig. 4, the SHLR-SV, which equipped with both the SPIRiT and virtual coil information, produces the image with lowest error (Fig. 4 (f)) while the basic SHLR and the SHLR with only one extra information enforcement yield results with obvious errors ( Fig. 4 (b-c) and (e)). The results indicate the effectiveness of utilizing the two information in parallel imaging. 3) Parameter-dimensional information: Extra information can be incorporated to improve reconstructions of multidimensional MRI, e.g., information from parameter imaging.
Conventional parameter imaging acquires a series of images using different parameters. The image intensity variation along the parameter dimension then is used in data fitting to estimate the tissue intrinsic parameters, such as longitudinal relaxation time (T1), transversal relaxation time (T2). The tissue parameters are important in diagnosing the diseases including nervous, musculoskeletal, liver, and myocardial [14], [30]- [34]. It is worthy to mention that the signal along the parameter dimension can be modeled as mono or few exponentials [31], [34], [35], which indicates the low-rankness of the signal [15], [35], [36]. The low-rank property along the parameter dimension can be utilized to enhance the reconstructions. Motivated by ALOHA [35], we reconstruct images in the frequency encoding dimension one-by-one (Fig. 5). The difference between ALOHA and our proposed method is that ALOHA constructed a block Hankel matrix from the data on the phase encoding-parameter (PE-P) plane while the proposed method converted rows or columns on the phase encodingparameter plane into Hankel matrices separately. The seperable low-rank Hankel matrix strategy reduces the reconstruction time to a certain extent (Fig. 6).
Let X Parameter ∈ C M ×N ×L×J denote the desired parameter imaging MRI image. Denote X PE-P m ∈ C N ×L×J as the image at the m-th position on the frequency encoding dimension of X Parameter . Let Y Parameter denote the acquired k-space data of parameter imaging, with zero-filling at unacquired position and 1D inverse Fourier transform on along the frequency encoding dimension. And let Y PE-P m represent the data at the m-th position on the frequency encoding dimension of Y Parameter . F PE denotes the 1D Fourier transform along the phase encoding dimension. The reconstruction model for parameter imaging is as follow To further reduce the reconstruction error, we establish a reconstruction approach that takes advantages of both virtual coil and parameter-dimension information. Notice that SPIRiT exploits the correlation of image of different coils, it appears not appropriate for the recovery of PE-P plane images. Thus, the proposed method stems from the basic model and introduces two prior information, smooth phase of the image and low-rankness in the parameter dimension: As shown in Fig. 6, the SHLR-VP which utlizes the virtual coil acheives the lowest error, and the computational time of the proposed SHLR-VP is reduced to 2/3 of that of ALOHA. The results demonstrate the advantages of the proposed method in terms of low reconstruction error and fast reconstruction speed for parameter imaging.

C. Numerical Algorithm
We adopted the alternating direction method of multiplier (ADMM) [37] to solve both the proposed SHLR-SV in Eq.
where ·, · represents the inner product of matrix, and β the penalty parameter. The solution of Eq. (9) can be obtained by alternatively solving the following sub-problems: where τ is the step size.
where the upper subscript * denotes the adjoint operator. Here, for any matrix A, vec −1 (vec (A)) = I, therefore, we omitted the writing of vec −1 (vec (·)) for simplicity in the following.
For fixed X (k+1) and D row where S 1/β (·) denotes the singular value threshold operator with the threshold of 1/β. Similarly, the Z col n (k+1) can also be obtained by The numerical algorithm for SHLR-SV is summarized in Supplementary Material S1.
2) SHLR-VP: Considering the space limitations, please refer to Supplementary Material S2 for detailed derivations for solving the SHLR-VP model.

V. RESULTS
In this section, we will evaluate the two proposed approaches, SHLR-SV and SHLR-VP using in vivo data from different MRI scanners.
In all experiments, the multi-coil images are combined by a square root of squares (SSOS). We adopted the RLNE [38] and mean structure similarity index measure (MSSIM) [39] as objective criteria to quantify the reconstruction performance. A lower RLNE, the higher consistency lies between the fully sampled image and the reconstruction image. A higher MSSIM indicates higher detail preservation in reconstruction. The calculation of RLNE and MSSIM can be found in Supplementary Material.

A. Parallel Imaging
In this subsection, we assessed the reconstruction results and running time of SHLR-SV and compared the results of SHLR-SV with three state-of-the-art reconstruction approaches including 1 -SPIRiT [25], AC-LORAKS [24], and STDLR-SPIRiT [22]. These three methods and SHLR-SV require ACS for reconstruction. More specifically, 1 -SPIRiT reconstructs image by exploiting kernel estimation and the sparsity of the image in the transform domain. Both AC-LORAKS and STDLR-SPIRiT utilize the low-rankness of block Hankel matrix (or structured matrix) of k-space to recover MRI image. In addition, AC-LORAKS utilizes ACS for reducing the algorithm complexity [24].
The codes of 1 -SPIRiT are shared online by Dr. Michael Lustig [40] and codes of AC-LORAKS [41] are shared at Dr. Justin P. Haldar's website. Here, we adopted the Sbased AC-LORAKS with virtual coil which makes use of phase constraints. The reason why we chose S-based is that LORAKS with phase constraint provides results with lower error compared to other constraints [18]. Parameters of all the compared methods are optimized to obtain the lowest RLNE. For the proposed methods, the effect of parameters setting will be discussed in Supplementary Material S5. Two brain datasets acquired from healthy volunteers are adopted in experiments. The dataset depicted in Fig. 7 (a) is obtained from a 3T SIEMENS MRI scanner (Siemens Healthcare, Erlangen, Germany) equipped with 32 coils using T2weighted turbo spin echo sequence (matrix size = 256 × 256,  Fig. 8 (a) is acquired from a 3T SIEMENS Trio whole-body scanner (Siemens Healthcare, Erlangen, Germany) equipped with 32 coils using the 2D T2-weighted turbo spin echo sequence (matrix size = 256×256, TR/TE = 6100/99 ms, FOV = 220 mm × 220 mm, slice thickness = 3 mm). Four virtual coils are compressed from the acquired data of 32 coils [42].
We first test the proposed method for parallel imaging using the widely-used sampling pattern adopted in commercial MRI scanners named uniform undersampling pattern. Both 1 -SPIRiT ( Fig. 7 (b)) and AC-LORAKS (Fig. 7 (c)) have strong undersampling artifacts. STDLR-SPIRiT ( Fig. 7 (d)) and SHLR-SV (Fig. 7 (e)) show the good ability of artifacts removing. But it is worthy to note that SHLR-SV provides the image with lower error than STDLR-SPIRiT.
We then test the proposed method using non-uniform undersampling patterns, including 1D Cartesian and 2D random undersampling patterns. In the reconstruction under Cartesian undersampling pattern, AC-LORAKS (Fig. 8 (c)) yields results exhibiting obvious artifacts inside the brain area, whereas the ringing artifacts also remain in the reconstructed image of 1 -SPIRiT ( Fig. 8 (b)). Both STDLR-SPIRiT (Fig. 8 (d)) and SHLR-SV (Fig. 8 (e)) provide the image with nice artifacts suppression, however, the reconstruction error of STDLR-SPIRiT appear slightly larger than that of SHLR-SV, especially in the skull area.
In addition, SHLR-SV allows promising reconstruction on the undersampling patterns with partial Fourier (PF), which is commonly adopted in commercial equipment for further acceleration [26], indicating the potential to be robust to patterns. As shown in Fig. 9, under the case with PF, the reconstruction errors of three compared approaches, including STDLR-SPIRiT, apparently increased, whereas SHLR-SV still  yields nice reconstruction with low reconstruction error. In the meanwhile, the proposed SHLR-SV outperforms other methods in terms of RLNE and MSSIM, indicating the excellent ability of artifact suppression and detail preservation. We include more data and sampling patterns to further evaluate the proposed SHLR-SV approach. The detailed results can be found in Supplementary Material S3. Similar phenomenons were observed in the results indicating SHLR-SV is robust to the different data from different subjects, and also robust to different sampling patterns. It can be seen that SHLR-SV outperforms the 1 -SPIRiT and AC-LORAKS, and shows comparable results with STDLR-SPIRiT. Given the patterns with PF, SHLR-SV gains a slight improvement than STDLR-SPIRiT with respect to RLNE and MSSIM. However, the computational time of the proposed method is significantly reduced compared to that of STDLR-SPIRiT (8× faster). Overall, the proposed method provides faithful reconstruction with acceptable computational time, having the potential to be used in clinical applications.

B. Parameter Imaging
In this subsection, we validate the proposed SHLR-VP method using a acquired in-vivo T2 mapping data. We compare the proposed method named SHLR-VP with two state-ofthe-art mapping reconstruction methods including MORASA [15], and ALOHA [35]. These two methods both utilized the property of exponential in the parameter dimension. MORASA modeled the time-domain signal into a exponential function and the enforced the sparsity of the image at each time point. ALOHA made use of the annihilation relationship of the data and lifted the signal into a higher dimensional space to ensure the low-rankness. The codes of MORASA are shared by Dr. Xi Peng and codes of ALOHA [43] are shared at Dr. Jong Chul Ye website.
A fully sampled brain dataset was acquired from a healthy volunteer on a 3T 12 coils MRI scanner (Siemens Healthcare, Erlangen, Germany) using turbo spin echo sequence (TR = 4000 ms, 15 TEs from 8.8 to 132 ms with 8.8 ms sapcing, FOV = 200 mm × 200 mm, matrix size = 192×192, slice thickness = 3 mm). To estimate the T2 map, standard nonlinear least square fitting method was performed in the selected region-ofinterest on a pixel-by-pixel basis. T2 relaxation values outside the reasonable range were excluded from the T2 map.
We performed reconstructions under two different reduction factors, R = 6 and R = 8. In both cases, the MORASA images appear apparent string-like artifacts (Fig. 10 (e) and (i)). Compared to ALOHA, SHLR-VP produces mapping results with better accuracy. As can be seen from the images, the mapping error is obviously smaller than that of ALOHA. The RLNE and MSSIM also indicate that SHLR-VP persevere the image structure better than the comparison methods.
In addition, we validated the proposed method using other datasets, and the sampling patterns with PF. Please refer to the Supplementary Material S4 for detailed results. The results demonstrate that the proposed method SHLR-VP yields reconstructions with the lowest error in both reconstructed images and the quantization T2 map.

A. The Number of ACS
The ACS is crucial for traditional parallel reconstruction methods. These methods estimate sensitivity maps or kernels from ACS, once the number of ACS lines is limited, the estimation turns inaccurate, leading to failures of reconstruction [19]. Four comparison methods for parallel imaging require ACS to reconstruct images, thus, in this subsection we discuss the effect of the number of ACS on the reconstructed results.
The results shown in Fig. 11 indicates the proposed SHLR-SV maintains the robustness to ACS lines. When ACS lines (> 14 lines) are relatively sufficient, all these comparison methods provide reconstruction with relatively low reconstruction error. Especially, the proposed method SHLR-SV achieved the lowest RLNE and the highest MSSIM. However, when the number of ACS lines is extremely low (= 8), 1 -SPIRiT and AC-LORAKS fail to produce satisfying results (the RLNE reaches closely to 0.2). In contrast, the STDLR-SPIRiT and SHLR-SV, still work well, even with only 8 ACS lines offering an RLNE of approximately 0.07. The observation demonstrates that the proposed SHLR-SV still inherits the robustness of structured low-rank methods.

B. Computational Time
All reconstruction software were ran on a CentOS 7 computation server with two Intel Xeon CPUs of 3.5 GHz and 112 GB RAM. The computational time is obtained by averaging the reconstruction time of 10 Monte Carlo tests. All comparison approaches were implemented in MATLAB.
As shown in Fig. 12, the running time of SHLR-SV is reduced to 1/8 of that of STDLR-SPIRiT, which substantially alleviates the burden of lengthy reconstruction time of STDLR-SPIRiT. Moreover, it can be observed that the proposed method runs relatively slower than 1 -SPIRiT and AC-LORAKS. However, the additional computational cost of SHLR-SV should be acceptable considering its improvement in image reconstruction.
The computational time on parameter imaging reconstruction is also depicted in Fig. 12. It shown that the proposed SHLR-VP reconstructs data at a faster speed than MORASA and ALOHA. Besides, it is worthy to notice that runtime of the proposed method can be further reduced with GPU or MEX/C implementations.

VII. CONCLUSION
In this work, we attempted to alleviate the huge computational complexity and lengthy reconstruction time problems of the state-of-the-art structural low-rank approaches. To reduce the computational time, we proposed a separable Hankel lowrank reconstruction method, named SHLR, to enforce the Hankel low-rankness of each row and each column of the signal of interest. But this way to enforce the low-rankness is sub-optimal compared to block Hankel methods, though holds the advantage of swift runtime. To enhance the performance of SHLR, we introduced prior information. For parallel imaging, we explored the correlation between rows/columns and the virtual coil technique, and proposed the SHLR-SV model. Regarding parameter imaging, we assumed that the signal intensity along the parameter dimension varies in an exponential way. We enforced the exponential property and virtual coil technique and proposed a model named SHLR-VP for accelerating parameter imaging.  This experiment is carried out with a series of Gaussian Cartesian sampling patterns but all at the same sampling rate of 0.34. The difference between each sampling pattern is the number of ACS lines. Note: all experiments are conducted on the brain data shown in Fig. 8 (a). All methods are performed with the same parameters respectively, with which these methods provide the lowest RLNE when ACS lines equal 22. Experimental in-vivo results showed that the proposed two approaches enable better results than the state-of-the-art methods. Notably, the proposed methods allow reconstructions with smaller errors and faster reconstruction speed. The numerical algorithm for SHLR-SV is summarized in Algorithm S1.
Update X (k+1) by solving Eq. (15); 3: Update Z row m (k+1) and Z col The solution of Eq. (S1) is obtained by alternatively solving the following sub-problems: The solution of Eq. (S2) Similarly, the Eq. (S4) can be solved by The numerical algorithm for SHLR-VP is summarized in Algorithm S2.

S3. MORE EXPERIMENTAL RESULTS FOR PARALLEL IMAGING
In this section, more undersampling patterns and brain images are adopted for further validation in parallel imaging. We adopted RLNE and MSSIM as evaluation criteria. The RLNE is calculated through: where x andx denote the vectorized fully sampled and reconstructed SSOS image, respectively. The mean measure of structural similarity (MSSIM) on two images A and B is calculated by where A and B represent fully sampled and reconstructed SSOS images, respectively; µ ai , µ bi , σ ai , σ bi and σ aibi respectively denote the means, standard deviations and covariance of the local window a i and b i ; R the number of local windows. Constants C 1 and C 2 are introduced to avoid the case when the denominator multiplied by µ 2 ai + µ 2 bi and σ 2 ai + σ 2 bi is close to zero. Here, we first validate the proposed method under 2D random undersampling pattern. The 2D random undersampling patterns, which have stronger incoherence, are adopted here to simulate the 3D imaging with two phase encoding. As shown in Figs. S1 (b-e), all these four approaches reconstruct with artifact-free images. The reconstruction image of AC-LORAKS seems relatively more noise, but the error distributes randomly inside the skull, being preferred in clinical application. Errors of 1 -SPIRiT, STDLR-SPIRiT and SHLR-VC have reached a promising low level, assuring reliable reconstructions. It can also be seen that the SHLR-VC and STLR-SPIRiT provides slightly better reconstruction than that of 1 -SPIRiT in terms of RLNE and MSSIM.
Similarly, SHLR-SV outperforms the comparison methods under 2D random undersampling with PF. As shown in Fig.  S2 (g-j), the proposed method yields the reconstruction with lowest error distribution. The metrics listed in Table S1 also demostrate the superiority of the proposed method.
Besides, we adopted two more in vivo brain data under six different undersampling patterns to demonstrate the performance of the proposed method (shown in Fig. S3). Two brain datasets acquired from healthy volunteer. The first dataset shown in Fig. S3 (a) is acquired from a 3T SIEMENS MRI scanner (Siemens Healthcare, Erlangen, Germany) equipped with 32 coils using T2-weighted turbo spin echo sequence (matrix size = 256 × 256, TR/TE = 6100/99 ms, FOV = 220 mm × 220 mm, slice thickness = 3 mm). Four virtual coils are compressed from the acquired data of 32 coils [1]. The other dataset shown in Fig. S3 (e) is obtained from a 3T SIEMENS MRI scanner (Siemens Healthcare, Erlangen, Germany) equipped with 32 coils using the 2D T1-weighted FLAIR sequence (matrix size = 256 × 256, TR/TE = 3900/9.3 ms, FOV = 200 mm × 200 mm, slice thickness = 5 mm). Eight virtual coils are compressed from the acquired data of 32 coils [1].
The quality metric of reconstruction results are listed in Table S1. The same phenomenon as that in the main text can be observed. The proposed method in parallel imaging, SHLR-SV, is superior to 1 -SPIRiT and AC-LORAKS, and comparable to the STDLR-SPIRiT, in terms of RLNE and MSSIM. When it turns to undersampling patterns with PF, SHLR-SV performs slightly better results than STDLR-SPIRiT does.

S4. MORE EXPERIMENTAL RESULTS FOR PARAMETER IMAGING
In this section, we adopted two in vivo T2 mapping data to verify the performance of the proposed method. We first test the comparison methods under the sampling pattern with partial Fourier. The dataset used here was fully sampled on a 3T SIEMENS MRI scanner (Siemens Healthcare, Erlangen, Germany) with 20 coil using turbo spin echo sequence (TR = 4000 ms, 15 TE with 8.2, 16,25,33,41,49,58,66,74,82,91,99,107,115, and 124 ms, FOV = 220 mm × 220 mm, matrix size = 192 × 192, slice thickness = 4 mm). We removed 4 coils of data with strong noise, and finally, 16 coils of data were used for reconstructions.
Look at the reconstructed image at 9-th echo, we can find that the MORASA image has obvious sampling artifacts ( Fig. S4  (b) and (e)). The artifacts can still be seen in its mapping results (Fig. S4 (i) and (I)). Both ALOHA and SHLR-VP remove the sampling artifacts promisingly (Fig. S4 (c) and (d)). However, we can still see significant difference between their error images ( Fig. S4 (f) and (g)). The SHLR-VP error is much smaller than ALOHA. The associated mapping result of SHLR is more consistent with that of ALOHA (Please see the zoom-in image in Fig. S4 (j) and (k)).
Besides, we conducted experiments on undersampling patterns with different acceleration factors and patterns with/without PF. Another slice of T2 mapping data shown in Fig. 10 in the main text was adopted here. The results are summerized in Table S2. As shown in Table S2, the proposed method for parameter imaging, SHLR-VP, yields reconstruction with lower error and better structural preservation than that of comparison methods in both reconstructed images and estimated T2 maps. Especially when the acceleration factor is higher (R = 8) or when the undersampling patterns with PF, the improvement gained by SHLR-VP appears more obvious.

S5. PARAMETER SELECTION OF THE PROPOSED APPROACHES
In this subsection, the effect of parameters setting of the two proposed methods will be discussed, including regularization parameters λ and λ 1 for parallel imaging and λ and λ 2 for parameter imaging. The brain image in Fig. 8 (a) and the sampling pattern in Fig. 8 (f) has been used to perform parallel imaging reconstructions, and the T2 mapping data shown in Fig. 10 with the R = 8 reduction factor has been adopted for parameter imaging reconstruction.
The reconstruction errors versus different values of regularization parameters of SHLR-SV are shown in Fig. S5 (a). As can be seen in Fig. S5 (a), there exists a wide range of λ 1 (10 0 ≤ λ 1 ≤ 10 4 ) and λ (10 3 ≤ λ ≤ 10 6 ) leading to relatively low reconstruction errors. Too small or too large values of λ and λ 2 produce higher reconstruction errors. Similar trend was observed in the parameter imaging results (Fig. S5 (b)) but the ranges of λ and λ 2 are smaller. There exists a range of λ 2 (1 ≤ λ 2 ≤ 5) and λ (10 3 ≤ λ ≤ 10 5 ) leading to relatively low reconstruction errors.