Optshrink LR + S: accelerated fMRI reconstruction using non-convex optimal singular value shrinkage

This paper presents a new accelerated fMRI reconstruction method, namely, OptShrink LR + S method that reconstructs undersampled fMRI data using a linear combination of low-rank and sparse components. The low-rank component has been estimated using non-convex optimal singular value shrinkage algorithm, while the sparse component has been estimated using convex l 1 minimization. The performance of the proposed method is compared with the existing state-of-the-art algorithms on real fMRI dataset. The proposed OptShrink LR + S method yields good qualitative and quantitative results.


Introduction
Functional magnetic resonance imaging (fMRI) is one of the most significant noninvasive and non-ionizing diagnostic imaging modality [1,2]. It measures blood oxygenated level dependent (BOLD) signal for localizing brain activity [3]. However, despite the advancements in fMRI scanners, one of the biggest limitations of fMRI modality is slow imaging compared to the other medical imaging modalities [4].
Conventionally, parallel imaging techniques such as Sensitivity Encoding (SENSE) [5,6], Generalized Autocalibrating Partially Parallel Acquisitions (GRAPPA) [7,8], and Simultaneous Acquisition of Spatial Harmonics (SMASH) [9] are used to accelerate magnetic resonance imaging (MRI). Here, the basic principle involves use of multiple receiver coils with complementary sensitivity information. SENSE, GRAPPA, or SMASH reconstruct MRI images from multiple k-space undersampled images acquired on different coils. For the case of fMRI, only kt GRAPPA is able to accurately reconstruct fMRI images [10]. However, this method introduces strong temporal autocorrelation in the data that limits the extent of undersampling of fMRI data [10].
Apart from parallel imaging, compressed sensing (CS)based fMRI reconstruction is another attractive method of fMRI acceleration [11][12][13][14][15][16][17]. Similar to parallel imaging, less data are acquired in the k-space (spatial frequency domain) in CS resulting in accelerated fMRI data acquisition. However, unlike GRAPPA and SENSE that does not exploit information contained across time frames, CS exploits information across time frames leading to sparse representation and hence provides good reconstruction quality. Reconstruction of full fMRI data from this less or undersampled data requires efficient reconstruction algorithm. Researchers have proposed various methods for efficient reconstruction from undersampled k-space measurements [11][12][13][14][15][16][17]. These methods largely rely on compressive sensing and reconstruct data using an optimization framework under certain constraints. Often, fMRI data are assumed to be sparse in some transform domain. Theoretical studies have shown that it is possible to recover sparse signals by l 1 norm minimization [18]. For example, in [13], undersampled fMRI data are reconstructed using CS with sparsity of fMRI data in the wavelet domain, wherein orthogonal Daubechies wavelet is used as the sparsifying basis. This is to note that CS-based sparse recovery methods are being used extensively in many applications including other medical imaging modalities [19,20] and in videos [21,22].
In general, fMRI data matrix X, i.e., one fMRI slice data stacked over time, is observed to be low rank. Hence, lowrank constraint can be imposed in the CS-based optimization framework to recover fMRI data slice by slice. Recently, k-t FASTER method has been proposed on similar lines that recovers fMRI signal via hard thresholding of singular values of low-rank data matrix X in the CS framework [11].
In [15], a new LR ? S method had been proposed to reconstruct fMRI data that uses a linear combination of low-rank and sparse components, i.e., This decomposition of data into low-rank and sparse component is popularly known as robust principal component analysis (RPCA) in the literature [23,24]. In RPCA, convex optimization-based approaches are used to recover low-rank and sparse components from matrix X. In [15], fMRI reconstruction is solved via iterative estimation of L and S components using convex optimization-based approaches. In noise-free scenarios, convex approaches may still provide reasonable solution for non-convex problems [25]. However, in noisy settings, such as in fMRI with low signal-to-noise ratio (SNR), convex optimizationbased methods may not provide optimum or close to optimum solution. For quality accelerated fMRI reconstruction in noisy settings, improved low-rank matrix and sparse matrix estimation are necessary. There has been a great interest to recover low-rank matrix from noisy measurements in various fields such as statistical signal processing [26][27][28], machine learning [29], and estimation and classification problems [30]. This motivates us to explore an improved method of accelerated fMRI reconstruction that can recover denoised low-rank matrix and sparse component from the undersampled k-space data.
We use optimal singular value shrinkage denoising algorithm (OptShrink), a data-driven method, recently used for denoising of low-rank matrix [31]. We call the proposed method as Optshrink LR ? S method. In [31], OptShrink has been shown to have improved performance over singular value thresholding (SVT) in the recovery of data with missing entries. The OptShrink method requires noisy low-rank matrix and its rank estimate as input and provides denoised low-rank matrix estimate.
The proposed Optshrink LR ? S fMRI reconstruction method is compared with other offline fMRI reconstruction methods such as direct inverse Fourier transform (IFT), LR ? S [15], and CS with wavelet sparsity [13] methods. We compare reconstruction results using different methods at both the subject-and group level at different acceleration factors. Our proposed OptShrink LR ? S method reconstructs fMRI data with greater accuracy compared to other methods even at lower sampling ratios.
The rest of this paper is organized as follows. Section 2 discusses fMRI reconstruction problem and presents the proposed Optshrink LR ? S reconstruction method. In Sect. 3, simulation results using Optshrink LR ? S and some of the existing methods are presented on real fMRI data. Conclusions are presented in the last section.

Materials and methods
In this section, we present the mathematical formulation of fMRI reconstruction problem followed by details of the proposed Optshrink LR ? S fMRI reconstruction method and description of the fMRI dataset used in simulations.

Problem formulation
The functional MRI imaging involves acquisition of contiguous brain slices over a number of time points. For each individual brain slice, Casorati matrix X 2 R nÂT is formed by stacking one brain slice over all time points [32], i.e., where T is the number of time points and n is the number of voxels in one brain slice. Hence, each column x i of X corresponds to data of a particular brain slice captured at one time point. Let us denote the undersampled k-space fMRI data of one brain slice captured over time by the matrix Y.
The relationship between undersampled k-space data Y and X is as follows: where F denotes the two-dimensional (2-D) Fourier transform operator, U is the measurement matrix that detects or captures fewer k-space measurements, and n 2 R nÂT denotes the measurement noise. In fMRI reconstruction problem, matrix X is required to be recovered from the undersampled k-space fMRI data measurements in matrix Y.

Reconstruction using low-rank plus sparse decomposition
In this paper, we are interested in accelerated fMRI data reconstruction using low-rank plus sparse decomposition.
Hence, in this subsection we first elucidate low-rank plus sparse reconstruction problem. Consider matrix X as the superposition of low-rank matrix L with rank m and sparse matrix S with sparsity s. Hence, matrix X can be denoted as X :¼ L þ S 2 R nÂT , where matrices L and S are required to be recovered, given a set of undersampled measurements Y and the corresponding measurement matrix U. The optimization problem for identifying matrixL 2 R nÂT and matrixŜ 2 R nÂT from Y and U can be written as: where the Frobenius norm k Y À UFðL þ SÞ k F is defined as k Y À UFðL þ SÞ k F ¼ Tr½ðY À UFðL þ SÞÞ T ðY À U FðL þ SÞÞ, TrðÁÞ denotes trace of matrix, and ðÁÞ T denotes transpose. Problem in (3) can be solved iteratively when there is incoherence between low-rank matrix L and sparse matrix S matrices [15,23,24]. It is observed that incoherence is guaranteed when low-rank matrix is not sparse and sparse matrix is not low rank [15,24].

Proposed Optshrink LR ? S method
In this subsection, we explain the proposed Optshrink LR ? S reconstruction method wherein the problem in (3) is solved by breaking it into two subproblems of estimating L and S as described below.

S subproblem
In (3), k S k 0 denotes l 0 norm that is equal to the number of nonzero values (= s) in matrix S. l 0 norm is a non-deterministic polynomial (NP) hard problem [33]. Thus, l 1 norm is generally used as the closest convex surrogate of l 0 norm [18,34]. l 1 norm is defined as absolute sum of values in matrix S and is used to obtain sparse solution [18,34]. Generally, soft thresholding (ST) is used to solve l 1 norm penalty on S defined as: where denotes component-wise product and k 1 is the soft thresholding regularization parameter on S. Recently in [15], sparse matrix is recovered using ST on S. We use similar approach of ST in this work to solve for sparsity on S.

L subproblem
Low-rank matrix recovery is ill-posed and NP hard [35]. One of the methods to solve this problem is via convex optimization using nuclear norm minimization [35].
Nuclear norm minimization implies l 1 penalty on singular values of matrix L that supports matrix L to be low rank. Global minimum of convex nuclear norm minimization is obtained by soft thresholding on singular values, known as singular value thresholding (SVT) [36].
To understand this, consider n Â T noisy low-rank matrix: where L is the noise-free low-rank matrix and d is a random noise matrix. Here, the goal is to estimate non-noisy low-rank matrix L from noisy matrixL. Let singular value decomposition (SVD) of matrixL is P q i¼1 r i u i v H i , where r i ; u i and v i are the singular values, left singular vectors, and right singular vectors, respectively; q ¼ minðn; TÞ denotes the rank ofL and ðÁÞ H denotes the conjugate transpose. Convex nuclear norm solution of (5) can recover non-noisy low-rank matrix via SVT [36] as: where definition of 'Soft' is same as defined in (4) and k 2 in (6) is the regularization parameter. Recently, in [15], low-rank matrix is recovered using SVT, where noisy input low-rank matrix is initialized from the previous iteration. The key idea behind SVT is to shrink nonsignificant singular values toward zero while keeping the singular vectors unchanged. However, nuclear norm minimization is an over-relaxing recovery solution of low-rank matrix [37].
In this paper, we propose to estimate non-noisy low-rank matrix or denoised approximation for the low-rank matrix in (5) that will provide an overall improved performance of fMRI signal reconstruction using low-rank plus sparse decomposition. In [31], best approximate noise-free lowrank matrix is obtained by optimal weighted combination of left and right singular vectors of input noisy matrixL in (5). Let us assume that low-rank matrix L has rank m [refer to (5)] and is given as where u i and v i are the left and right singular vectors of noisy matrixL and w i are unknown singular values. In order to recover L from the noisy matrixL in (5), the problem is formulated as: Using (7), we can rewrite (8) as: Optshrink LR ? S: accelerated fMRI reconstruction using non-convex optimal singular value... 67 where, l 0 norm on w in (9) signifies number of nonzero singular values equal to the rank m. In the above equation, singular vectors u i and v i are known and estimated using SVD of input matrixL. The closed form solution of singular values in (9) for every 1 i m is expressed as [31]: where R 2 R nÂT and is equal to diagðr mþ1 ; . . .r q Þ, q ¼ minðn; TÞ, r i denotes the ith singular value ofL, TrðÁÞ is equal to the trace of a matrix, and I is the identity matrix.
DðÁÞ is the D-transform which is defined as and D 0 ðÁÞ is defined as This algorithm is named as Optshrink [31]. OptShrink is a data-driven method, recently used for denoising of lowrank matrix in an application of signal recovery in missing data. It considers noisy low-rank matrix and its rank estimate [= m in (7)] as input, and provides denoised estimate of the low-rank matrix as the output. It is a non-convex solution that does weighing of singular vectors. It shrinks the corresponding singular values using truncated singular value decomposition (TSVD) and hence is called nonconvex optimal SVT. This algorithm works better than SVT [31].
Also, in [31], it has been shown that the solution of Optshrink is quite robust to input rank specification, and hence a rough estimate of rank [= m in (7)] at the input is sufficient. Another advantage of Optshrink is that there is no need to specify shrinkage parameter as is required in SVT [refer to k 2 in (6)]. In SVT, we need to tune k 2 for every dataset. It has been observed that Optshrink always outperforms SVT in the estimation of low-rank matrix.
In this paper, we propose to apply OptShrink for lowrank matrix estimation in fMRI reconstruction using lowrank plus sparse decomposition. fMRI data inherently have low signal-to-noise ratio (SNR). Hence, fMRI reconstruction with OptShrink for denoised low-rank matrix estimation should outperform existing low-rank plus sparse fMRI reconstruction method [15].

Overall solution of (3) using Optshrink LR ? S
Finally, Eq. (3) is solved iteratively using the proposed Optshrink LR ? S method as below: where A ¼ UF in (13) and j denotes an iteration number.
Here, ST is used to recover sparse matrix S as explained in Sect. 2.3.1 and Optshrink algorithm is used to solve for low-rank matrix L as explained in Sect. 2.3.2. Voxel time series are observed to be sparse in the Fourier domain. Hence, variation along rows of matrix X is assumed to be sparse in the Fourier domain. W in Algorithm 1 is the sparsifying matrix for the Fourier domain where Fourier transform is to be taken along the rows. Solution is updated at each iteration j. Algorithm 1 presents the pseudo-code of Optshrink LR ? S method.
Please note that low-rank component represents the background information that is highly correlated across data captured at different time points and sparse component represents the dynamic and uncorrelated counterpart.

Dataset description
To assess the performance of reconstruction methods, we have used two fMRI datasets in this paper: (1) Task-based fMRI dataset with false belief task (OpenfMRI publicly available dataset) 1 and (2) Resting-state Baltimore fMRI dataset (1000 functional Connectomes Project Data). 2

Task-based dataset
This dataset consists of acquisition of 36 axial interleaved brain slices with dimensions 72x72 at each time point with echo time (TE) equal to 35 ms and repetition time (TR) equal to 2 s [38]. These data are collected over 179 time points, resulting in the matrix X of size 5184 Â 179 for one brain slice. During the false belief experiment, the subject had to answer questions about stories that referred either to a person's false belief (mental trials) or to outdated physical representations such as an old photograph. For more details on this dataset, please refer to [38].

Resting-state dataset
These data are publicly available as part of the 1000 Functional Connectomes Project. This is a collection of resting-state fMRI dataset from a number of laboratories around the world. We use Baltimore resting fMRI data. This dataset consists of 23 subjects resting-state fMRI data, aged between 20 and 40 years of age, acquired while subjects' eyes were open and fixated on a screen. The repetition time (TR) is 2.5 s, size of a brain volume at one time point is 96 Â 96 Â 47, and the total no. of time points over which data are captured is 123.

Simulation results
Since both resting-state fMRI dataset (Baltimore dataset) and task-based fMRI dataset (false belief fMRI dataset) are fully sampled, we simulated undersampled k-space data Y in (2) by computing the Fourier transform of Casorati matrix X and then, retrospectively undersampling in the k-space using measurement matrix U. This measurement matrix is generated using radial sampling patterns. We used three radial sampling patterns with different acceleration factors for testing reconstruction performance: 6 radial lines, 12 radial lines, and 24 radial lines as described in [39]. Radial sampling pattern is chosen because this is one of the fastest kspace sampling methods in real-time application [39]. These radial sampling patterns sample more data in the lowfrequency region compared to the high-frequency region. This is to note that we have illustrated undersampling of fMRI data by retrospective sampling on the Cartesian grid because it allows sampling patterns to maintain incoherency among the columns of matrix X [39][40][41]. However, radial-Cartesian sampling grid is more realistic from the point of view of actual data acquisition [10,42]. Similarly in [12], variable density spiral sampling pattern has been used in MRI scanner and is shown to be robust against motion, off resonance, and gradients artifacts in compressed sensing fMRI application. However, our work is focused on development of robust reconstruction algorithm. This is to note that the proposed Optshrink LR ? S method reconstructs fMRI data as a superposition of low-rank and sparse matrix, where the lowrank component represents background information that is highly correlated across data captured at different time points and sparse component represents the dynamic and uncorrelated part. Since these assumptions are characteristic of fMRI data, they will remain valid irrespective of the sampling strategy used. Hence, although the proposed work is general and can be used with any sampling pattern provided sampling incoherence is maintained, we project the use of realistic sampling patterns as the future work.
Data obtained from the database are called as original data in the manuscript. k-space data are acquired by considering 2-D Fourier transform of this original data. Since the original data are real and are provided without any phase information, we only considered the magnitude part of the reconstructed data. Thus, no assumption is made about the phase part of the data. This is to clarify that this is a standard method of testing newer MRI/fMRI reconstruction algorithms via simulation results.

Comparison with different methods
In this section, we provide results on fMRI reconstruction from undersampled k-space fMRI data using the proposed Optshrink LR ? S method, existing LR ? S method [15], direct inverse Fourier transform-based reconstruction, and reconstruction using CS with wavelet sparsity [13]. Below, 1 https://openfmri.org/dataset/. 2 http://www.nitrc.org/frs/?group_id=296 Optshrink LR ? S: accelerated fMRI reconstruction using non-convex optimal singular value... 69 we present a brief overview of each of these existing reconstruction methods used.
3.1.1 Low-rank plus sparse (LR ? S) method [15] This method reconstructs fMRI data using superposition of low-rank and sparse matrix components [15], and hence, the optimization problem is:L We empirically selected k 1 ¼ 2 and k 2 ¼ 200 in the above Eq. (14) using the L-curve method [43]. Minimum normalized mean square error (NMSE) is obtained in the L-curve at the above chosen ks for the existing LR ? S method. This is to note that we used same values of ks in the proposed Optshrink LR ? S method. Thus, the values of ks are optimally selected for the existing LR ? S method and not for the proposed Optshrink method for presenting the comparative results.

Direct IFT
This method computes 2-D inverse Fourier transform (IFT) of given k-space fMRI data Y and reconstructs X as shown below:

CS with wavelet sparsity (CSWD) [13]
Wavelet sparsity assumes signal to be sparse in the wavelet domain [13], and hence, the optimization problem is:  Optshrink LR ? S: accelerated fMRI reconstruction using non-convex optimal singular value... 71 where W is a wavelet matrix operator. We used Daubechies' orthogonal wavelet 'db4' (filter lengths 8) with three-level decomposition as the sparsifying basis to exploit wavelet sparsity as used in [13]. This method requires one parameter k 3 to be specified as shown in (16).
In [44], k 3 is restricted to satisfy the below condition: In order to meet the above condition, we chose that meets (17). For all reconstruction methods, we set the maximum number of iterations equal to 500 and use the following convergence criterion: objective function valueðendÞÀ objective function valueðend À 1Þ\10 À5 , also specified in Algorithm 1. Figure 2 shows objective function value versus number of iterations on one subject of the taskbased dataset with the proposed Optshrink LR ? S method. From this figure, we observe that the objective function converges monotonically. We observed the same  trend with every data, and hence, we may safely state that the proposed Optshrink LR ? S method convergences to provide solution. Figure 3 shows one of the reconstructed brain slices of task-based fMRI dataset (false belief task) using the proposed Optshrink LR ? S and the existing LR ? S [15] method. Reconstructed data are visually shown at different radial lines in Fig. 3a, b, c corresponding to the middle slice (slice number 18 of 36 no. of total slices) captured at the 100th time point.
Likewise, Fig. 4 shows one of the reconstructed brain slices of resting-state dataset (Baltimore dataset) using the proposed Optshrink LR ? S and the existing LR ? S [15] method. Reconstructed data are visually shown at different radial lines in Fig. 4a, b, c corresponding to the middle slice (slice number 24 of 47 no. of total slices) captured at the 100th time point.
Following observations can be drawn from the reconstructed slices of both task-based and resting-state data shown in Figs  Task-based data-false belief task fMRI data, subject no. 1 Fully sampled fMRI data-cluster size = 26, Z score All the above observations indicate that we can reconstruct fMRI data with greater quality by sampling much lesser measurements in k-t space with the proposed Optshrink LR ? S method compared to the existing LR ? S method. Hence, higher acceleration rate is possible with Optshrink LR ? S method. Figures 5 and 6 show quantitative results via normalized mean square error (NMSE) versus time for both the dataset with different radial sampling patterns where: jj Á jj 2 denotes l 2 norm and, I andÎ are the original and reconstructed brain slices, respectively. NMSE values are computed between reconstructed and original slice at each time point. We represent reconstructed results using LR ? S method and Optshrink LR ? S method. In consonance with the qualitative results of Figs. 3 and 4, we observe higher NMSE with LR ? S method compared to the proposed Optshrink LR ? S method. While the NMSE increases rapidly with decrease in radial lines with LR ? S method, it remains quite consistent with Optshrink LR ? S method. In order to assess other reconstruction methods quantitatively, we present reconstruction results in Table 1 in terms of NMSE and peak signal-to-noise ratio (PSNR) on both the datasets. From Table 1, we note that NMSE increases with decrease in the number of radial lines, i.e., with fewer k-space measurements with existing reconstruction methods. On the other hand, the proposed Optshrink LR ? S method shows consistent PSNR and NMSE values at different radial lines.   Fig. 12 False belief fMRI data shown on sagittal, coronal, and axial planes: a fully sampled fMRI data; b smoothed fully sampled fMRI data; c reconstructed fMRI data using LR ? S (6 radial lines) (without smoothing); d reconstructed fMRI data using LR ? S (6 radial lines) (with smoothing); e reconstructed fMRI data using proposed Optshrink LR ? S method (rank = 1) (6 radial lines) (without smoothing); f reconstructed fMRI data using proposed Optshrink LR ? S method (rank = 1) (6 radial lines) (with smoothing) Moreover, Optshrink LR ? S results are similar with different input rank specification. These results are in line with the qualitative results observed with the reconstructed slice quality in Figs. 3 and 4. This is to note that the proposed Optshrink LR ? S method estimates a denoised version of low-rank matrix and hence yields better results.
In order to further ascertain the robustness of Optshrink LR ? S method to rank, we present NMSE versus rank and PSNR versus rank for both the dataset in Figs. 7 and 8. We use 6 radial lines for undersampling k-space data and provide different rank as input to Optshrink LR ? S method. The reconstruction accuracy remains similar for different rank values, and hence, any rough estimate of rank may be provided as input to this method for fMRI signal reconstruction.

Group-level analysis
In Figs. 9 and 10, we present NMSE and PSNR results for five subjects each of task-based dataset and resting-state dataset, respectively. We use 6 radial lines for undersampling the k-space data. These results indicate that our proposed Optshrink LR ? S is robust to subject variability and are reproducible across subjects.

Subject-level statistical analysis on activation maps
In this section, we would like to study the effectiveness of Optshrink LR ? S method with reference to brain activation detection. To this end, reconstruction is performed on the task-based dataset (false belief task) using LR ? S method and Optshrink LR ? S (rank = 1) method. Preprocessing of the original and the reconstructed fMRI dataset are performed using SPM12. 3 We performed motion correction that is used to suppress motion-related artifacts. In general, motion correction is followed by smoothing as a preprocessing step so that the noise is Gaussian-distributed (by Central Limit Theorem). This establishes the validity of statistical tests using general linear model (GLM)-based analysis, a univariate approach, used for detecting brain activation in task-based fMRI data [45]. Since Optshrink LR ? S method is supposed to provide denoised low-rank matrix, we tested the robustness of the proposed method on brain activation detection both with and without smoothing in the preprocessing pipeline.
In GLM, linear model is fitted to each voxel time series using the design matrix corresponding to the task stimulus. The estimated parameters are used to build statistical parametric maps (SPMs) [46]. Figure 11 shows the design matrix for the false belief dataset that consists of five conditions. First four conditions correspond to four different block stimuli (false belief story, false belief question, false photograph story, false photograph question [38]) that are convolved with the canonical hemodynamic response function (HRF) and form first four columns of the design matrix. The last column captures the linear trend of data.
Reconstruction is performed on undersampled fMRI data on three radial sampling patterns of 6, 12, and 24 radial lines. Figures 12, 13, and 14 show the corresponding statistical maps obtained using (a) original fully sampled k-t space data without smoothing operation (b) original smoothed fully sampled data, (c) reconstructed data using LR ? S method (without smoothing), (d) reconstructed data using LR ? S method (with smoothing), (e) reconstructed data using the proposed Optshrink LR ? S method (without smoothing), and (f) reconstructed data using the proposed Optshrink LR ? S method (with smoothing). We present results on representative slices having peak voxel of activation, whereas Montreal Neurological Institute (MNI) position of this most active voxel is listed in Table 2. We also report cluster sizes and maximum z-scores values in this table. Activation maps are thresholded t test at cluster level with uncorrected p value = 0.05. Clusters with less than 12 voxels are rejected.
Brain activation maps using original fully sampled smoothed data show better results compared to the original fully sampled data without smoothing [compare (b) with (a) in Figs. 12,13,14]. As evident from Figs. 12, 13, 14 and Table 2, we notice that LR ? S method provides inferior results, while activation maps using data reconstructed with Optshrink LR ? S method (with smoothing in preprocessing) provides activation maps similar to those of (b). The MNI position of the most active voxel on the reconstructed data using the proposed Optshrink LR ? S method (with smoothing) is same as that obtained with the original data. Smoothing helps in increasing the sensitivity of BOLD time series. It can be noticed via increase in cluster size containing active voxels. Original smoothed fMRI data cluster size is 112 while without smoothed cluster size is 26. In the case of the proposed Optshrink b Fig. 13 False belief fMRI data shown on sagittal, coronal, and axial planes: a fully sampled fMRI data; b smoothed fully sampled fMRI data; c reconstructed fMRI data using LR ? S (12 radial lines) (without smoothing); d reconstructed fMRI data using LR ? S (12 radial lines) (with smoothing); e reconstructed fMRI data using proposed Optshrink LR ? S method (rank = 1) (12 radial lines) (without smoothing); (f) reconstructed fMRI data using proposed Optshrink LR ? S method (rank = 1) (12 radial lines) (with smoothing) LR ? S method (with smoothing), the cluster size of reconstructed data is 204. Also, these clusters of activation maps are consistently good at all acceleration factors. This clearly shows that the proposed method provides enhanced brain activation maps and is indeed better compared to other reconstruction methods.

Reproducibility of resting-state networks
In this section, we test the efficacy of the proposed Optshrink LR ? S method on resting-state fMRI dataset. We evaluate performance in terms of reproducibility of restingstate networks (RSNs). We compare RSNs of Optshrink LR ? S-based reconstructed data with the RSNs obtained from the fully sampled original fMRI data that is considered as the ground truth. RSNs are identified using the spatial independent component analysis (ICA) of GIFT toolbox. 4 Before ICA is applied, data are preprocessed. The first five fMRI brain volumes are discarded followed by slicetime correction. Next, realignment is done for motion correction followed by spatial normalization onto the Montreal Neurological Institute (MNI) space (3-mm isotropic voxels). In the end, brain volumes are spatially smoothed with a Gaussian kernel [full width half maximum (FWHM) = 6 mm].
After preprocessing, we utilize InfomaxICA algorithm in GIFT to obtain 100 independent spatial components. We identified 54 RSNs from mean maps of all five fully sampled ground truth fMRI data (corresponding to five subjects) after removing artifact components. These RSNs can be broadly categorized into 10 RSNs: (1) Visual network (VN), (2) Somatomotor network (SMN), (3) Limbic network (LN), (4) Dorsal attention network (DAN), (5) Ventral attention network (VAN), (6) Default mode network (DMN), (7) Frontoparietal network (FPN), (8) Temporal ? Frontal network (TFN), (9) Subcortical network (SCN), and (10) Cerebellar network (CN). We also ran ICA on the reconstructed data. These dataset are reconstructed using 16.49% (12 radial lines) acquired samples in k-space using Optshrink LR ? S with rank one. We identified 56 RSNs from mean spatial components. These RSNs can be further classified into various categories as mentioned above. Figures 15 and 16 show some of the RSNs obtained from fully sampled ground truth data and Optshrink LR ? S reconstructed data. From this figure, we observe that RSNs identified by Optshrink LR ? S reconstructed data resemble with the ground truth fully sampled data. This shows that Optshrink LR ? S method is able to preserve functional characteristics of data. This is the most desirable need in neuroimaging research. Please note that we also ran ICA on reconstructed data using LR ? S. We observed more artifact components with a total 100 spatial components. This again validates our claim that the proposed Optshrink LR ? S method is working better than existing methods in the literature.

Conclusion
In this paper, we have proposed a new accelerated fMRI method, named Optshrink LR ? S method, for fMRI reconstruction from undersampled k-t space data. The proposed method exploits sparsity and low-rank decomposition with denoising to improve fMRI reconstruction accuracy. Comparison results demonstrate that the reconstruction performance of the proposed Optshrink LR ? S method is superior to existing methods at various acceleration factors. While the performance of the existing methods falls rapidly at faster acceleration rates, Optshrink LR ? S method performs consistently. Quantitative and qualitative results, group-level and subjectlevel analyses, show the superior performance of the proposed method. In addition, Optshrink LR ? S method provides enhanced brain activation maps that is an added but most useful advantage of the proposed method. MATLAB implementation of proposed algorithm is available online. 5 4 https://www.nitrc.org/projects/gift. 5 http://in.mathworks.com/matlabcentral/fileexchange/60836-optshrinklr?s-accelerated-fmri-reconstruction-using-non-convex-optimal-singularvalue-shrinkage. b Fig. 14 False belief fMRI data shown on sagittal, coronal, and axial planes: a fully sampled fMRI data; b smoothed fully sampled fMRI data; c reconstructed fMRI data using LR ? S (24 radial lines) (without smoothing); d reconstructed fMRI data using LR ? S (24 radial lines) (with smoothing); e reconstructed fMRI data using proposed Optshrink LR ? S method (rank = 1) (24 radial lines) (without smoothing); f reconstructed fMRI data using proposed Optshrink LR ? S method (rank = 1) (24 radial lines) (with smoothing) Fig. 15 Axial view of spatial maps of various RSNs where left part of each figure is from the original fully available dataset and right part is from the Optshrink LR ? S reconstructed data