Search tips
Search criteria 


Logo of biBrain Informatics
Brain Inform. 2017 March; 4(1): 65–83.
Published online 2017 January 10. doi:  10.1007/s40708-016-0059-x
PMCID: PMC5319953

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.

Keywords: Accelerated functional MRI, Low-rank recovery, Sparse recovery, Compressed sensing, k–t acceleration, Undersampling


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 [1117]. 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 [1117]. 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, low-rank constraint can be imposed in the CS-based optimization framework to recover fMRI data slice by slice. Recently, kt 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.,

XL(low rank) + S(sparse).

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 optimization-based 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 [2628], 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 ∈ ℝn×T is formed by stacking one brain slice over all time points [32], i.e., X = {xii = 1, …, T}, where T is the number of time points and n is the number of voxels in one brain slice. Hence, each column xi 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:

Y = ΦFXξ

where F denotes the two-dimensional (2-D) Fourier transform operator, Φ is the measurement matrix that detects or captures fewer k-space measurements, and ξ ∈ ℝ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: = LS ∈ ℝn×T, where matrices L and S are required to be recovered, given a set of undersampled measurements Y and the corresponding measurement matrix Φ. The optimization problem for identifying matrix L^Rn×T and matrix S^Rn×T from Y and Φ can be written as:

L^,S^=arg minL,S:rank(L)m,S0sY-ΦF(L+S)F,

where the Frobenius norm Y - ΦF(LS)‖F is defined as Y - ΦF(LS)‖F = Tr[(YF(L+S))T(Y - ΦF(LS))], 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), S0 denotes l0 norm that is equal to the number of nonzero values (= s) in matrix S. l0 norm is a non-deterministic polynomial (NP) hard problem [33]. Thus, l1 norm is generally used as the closest convex surrogate of l0 norm [18, 34]. l1 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 l1 norm penalty on S defined as:


where  ⊗  denotes component-wise product and λ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 l1 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 δ is a random noise matrix. Here, the goal is to estimate non-noisy low-rank matrix L from noisy matrix L~.

Let singular value decomposition (SVD) of matrix L~ is i=1qσiuiviH, where σiui and vi are the singular values, left singular vectors, and right singular vectors, respectively; q = min(nT) denotes the rank of L~ 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 λ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 low-rank matrix is obtained by optimal weighted combination of left and right singular vectors of input noisy matrix L~ in (5). Let us assume that low-rank matrix L has rank m [refer to (5)] and is given as


where ui and vi are the left and right singular vectors of noisy matrix L~ and wi are unknown singular values. In order to recover L from the noisy matrix L~ in (5), the problem is formulated as:

L^=arg minLL~-LFwithrank(L)=m.

Using (7), we can rewrite (8) as:

wopt=arg minwl0=mL~-i=1mwiuiviHF.

where, l0 norm on w in (9) signifies number of nonzero singular values equal to the rank m. In the above equation, singular vectors ui and vi are known and estimated using SVD of input matrix L~. The closed form solution of singular values in (9) for every 1 ≤ i ≤ m is expressed as [31]:


where Σ ∈ ℝn×T and is equal to diag(σm+1, …σq), q = min (nT), σi denotes the ith singular value of L~, 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( · ) is defined as


This algorithm is named as Optshrink [31]. OptShrink is a data-driven method, recently used for denoising of low-rank 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 non-convex 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 λ2 in (6)]. In SVT, we need to tune λ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 low-rank matrix estimation in fMRI reconstruction using low-rank 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 = ΦF 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. Ψ 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.

An external file that holds a picture, illustration, etc.
Object name is 40708_2016_59_Figa_HTML.jpg

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 Φ. 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 k-space sampling methods in real-time application [39]. Figure 1 shows these radial sampling measurement patterns. These radial sampling patterns sample more data in the low-frequency region compared to the high-frequency region.

Fig. 1
Radial sampling pattern on one slice: a 6 radial lines (12.856 acceleration factor); b 12 radial lines (6.065 acceleration factor); c 24 radial lines (3.495 acceleration factor)

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 [3941]. 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 low-rank 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, we present a brief overview of each of these existing reconstruction methods used.

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^,S^=arg minL,SY-ΦF(L+S)F2+λ1ΨS1+λ2L.

We empirically selected λ1 = 2 and λ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 λs for the existing LR + S method. This is to note that we used same values of λs in the proposed Optshrink LR + S method. Thus, the values of λs 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:

X^=arg minXY-ΦFXF2+λ3WX1,

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 λ3 to be specified as shown in (16). In [44], λ3 is restricted to satisfy the below condition:

λ3 < max (ΦT(IFT(Y))).

In order to meet the above condition, we chose

λ3 = 0.009 × max (ΦT(IFT(Y)))

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 task-based 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.

Fig. 2
Objective function value versus number of iterations

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.

Fig. 3
Task-based fMRI data–original and reconstructed slice no. 18 [left to right Original; LR + S; Optshrink LR + S (rank = 1); Optshrink LR + S (rank = 2); Optshrink LR + S ...

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.

Fig. 4
Resting-state fMRI data—original and reconstructed slice no. 24 [left to right Original; LR + S; Optshrink LR + S (rank = 1); Optshrink LR + S (rank = 2); Optshrink ...

Following observations can be drawn from the reconstructed slices of both task-based and resting-state data shown in Figs. 3 and 4:

  1. Slices reconstructed using LR + S method show a decline in quality with decrease in the number of radial sampling lines. On the other hand, reconstruction results with the proposed Optshrink LR + S method are quite consistent and the reconstruction quality does not fall by a great deal with the reduction in number of sampling lines.
  2. Slices reconstructed using different radial sampling patterns consistently show that LR + S method output is blurred and the slices have artifacts at the center and the boundary compared to the proposed Optshrink LR + S. This observation indicates that there is SNR loss with LR + S method that may lead to incorrect brain activation detection. On the other hand, slices reconstructed using the proposed Optshrink LR + S method are very clear and free of artifacts.
  3. Reconstruction using Optshrink LR + S method is robust to input rank specification. Hence, rank definition is not a bottleneck in the proposed Optshrink LR + S.

All the above observations indicate that we can reconstruct fMRI data with greater quality by sampling much lesser measurements in kt 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 and66 show quantitative results via normalized mean square error (NMSE) versus time for both the dataset with different radial sampling patterns where:


|| · ||2 denotes l2 norm and, I and I^ 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. 5
Normalized mean square error versus time points on task-based fMRI dataset (slice no. 18): a 6 radial lines (12.856 acceleration factor); b 12 radial lines (6.065 acceleration factor); c 24 radial lines (3.495 acceleration factor)
Fig. 6
Normalized mean square error versus time points on resting-state fMRI dataset (slice no. 24): a 6 radial lines (12.856 acceleration factor); b 12 radial lines (6.065 acceleration factor); c 24 radial lines (3.495 acceleration factor)
Table 1
Reconstruction results with different methods on two datasets

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.

Fig. 7
Task-based fMRI dataset, slice no. 18, time point 100: NMSE versus rank of the proposed Optshrink LR + S method using 6 radial lines
Fig. 8
Resting-state fMRI dataset—slice no. 24, time point 100: NMSE versus rank of the proposed Optshrink LR + S method using 6 radial lines

Group-level analysis

In Figs. 9 and and10,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.

Fig. 9
Task-based fMRI dataset—slice no. 18: NMSE versus subject number with the proposed Optshrink LR + S method (rank = 1) using 6 radial lines
Fig. 10
Resting-state fMRI dataset—slice no. 24: NMSE versus subject number with the proposed Optshrink LR + S method (rank = 1) using 6 radial lines

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.

Fig. 11
Design matrix of task-based fMRI dataset (false belief task)

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 kt 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.

Table 2
Statistical analysis results for uncorrected
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 ...
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 ...
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 ...

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,13, ,14].14]. As evident from Figs. 12, ,13,13, ,1414 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 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 resting-state 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 slice-time 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 and1616 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.

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
Fig. 16
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


In this paper, we have proposed a new accelerated fMRI method, named Optshrink LR + S method, for fMRI reconstruction from undersampled kt 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 subject-level 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


The first author would like to thank Visvesvaraya research fellowship, Department of Electronics and Information Technology, Ministry of Communication and IT, Government of India, for providing financial support for this work.

Compliance with ethical standards

Compliance with ethical standards

Conflicts of interest

On behalf of all authors, the corresponding author states that there is no conflict of interest.

Contributor Information

Priya Aggarwal,

Parth Shrivastava,

Tanay Kabra,

Anubha Gupta,


1. Feinberg DA, Yacoub E. The rapid development of high speed, resolution and precision in fMRI. NeuroImage. 2012;62(2):720–725. doi: 10.1016/j.neuroimage.2012.01.049. [PMC free article] [PubMed] [Cross Ref]
2. Ogawa S, Lee TM, Kay AR, Tank DW. Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proc Natl Acad Sci USA. 1990;87(24):9868–9872. doi: 10.1073/pnas.87.24.9868. [PubMed] [Cross Ref]
3. Worsley K, Friston K. Analysis of fMRI time-series revisited—again. NeuroImage. 1995;2(3):173–181. doi: 10.1006/nimg.1995.1023. [PubMed] [Cross Ref]
4. Frank LR, Buxton RB, Wong EC. Estimation of respiration-induced noise fluctuations from undersampled multislice fMRI data. Magn Reson Med. 2001;45(4):635–644. doi: 10.1002/mrm.1086. [PubMed] [Cross Ref]
5. Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P, et al. Sense: sensitivity encoding for fast MRI. Magn Reson Med. 1999;42(5):952–962. doi: 10.1002/(SICI)1522-2594(199911)42:5<952::AID-MRM16>3.0.CO;2-S. [PubMed] [Cross Ref]
6. Tsao J, Boesiger P, Pruessmann KP. kt blast and kt sense: dynamic MRI with high frame rate exploiting spatiotemporal correlations. Magn Reson Med. 2003;50(5):1031–1042. doi: 10.1002/mrm.10611. [PubMed] [Cross Ref]
7. Griswold MA, Jakob PM, Heidemann RM, Nittka M, Jellus V, Wang J, Kiefer B, Haase A. Generalized autocalibrating partially parallel acquisitions (grappa) Magn Reson Med. 2002;47(6):1202–1210. doi: 10.1002/mrm.10171. [PubMed] [Cross Ref]
8. Huang F, Akao J, Vijayakumar S, Duensing GR, Limkeman M. kt grappa: a k-space implementation for dynamic MRI with high reduction factor. Magn Reson Med. 2005;54(5):1172–1184. doi: 10.1002/mrm.20641. [PubMed] [Cross Ref]
9. Sodickson DK, Manning WJ. Simultaneous acquisition of spatial harmonics (smash): fast imaging with radiofrequency coil arrays. Magn Reson Med. 1997;38(4):591–603. doi: 10.1002/mrm.1910380414. [PubMed] [Cross Ref]
10. Chiew M, Graedel NN, McNab JA, Smith SM, Miller KL. Accelerating functional MRI using fixed-rank approximations and radial-cartesian sampling. Magn Reson Med. 2016;76:1825–1836. doi: 10.1002/mrm.26079. [PMC free article] [PubMed] [Cross Ref]
11. Chiew M, Smith SM, Koopmans PJ, Graedel NN, Blumensath T, Miller KL. kt faster: acceleration of functional MRI data acquisition using low rank constraints. Magn Reson Med. 2015;74(2):353–364. doi: 10.1002/mrm.25395. [PMC free article] [PubMed] [Cross Ref]
12. Fang Z, Van Le N, Choy M, Lee JH. High spatial resolution compressed sensing (hsparse) functional MRI. Magn Reson Med. 2015;76:440–455. doi: 10.1002/mrm.25854. [PMC free article] [PubMed] [Cross Ref]
13. Holland DJ, Liu C, Song X, Mazerolle EL, Stevens MT, Sederman AJ, Gladden LF, D’Arcy RCN, Bowen CV, Beyea SD. Compressed sensing reconstruction improves sensitivity of variable density spiral fMRI. Magn Reson Med. 2013;70(6):1634–1643. doi: 10.1002/mrm.24621. [PubMed] [Cross Ref]
14. Lu W, Li T, Atkinson I, Vaswani N (2011) Modified-cs-residual for recursive reconstruction of highly undersampled functional MRI sequences. In: Image Processing (ICIP), 2011 18th IEEE International Conference on, pp 2689–2692. doi:10.1109/ICIP.2011.6116222
15. Singh V, Tewfik A, Ress D (2015) Under-sampled functional MRI using low-rank plus sparse matrix decomposition. In: Acoustics, Speech and Signal Processing (ICASSP), 2015 IEEE International Conference on, pp 897–901. doi:10.1109/ICASSP.2015.7178099
16. Yan S, Nie L, Wu C, Guo Y. Linear dynamic sparse modelling for functional MR imaging. Brain Inform. 2014;1(1–4):11–18. doi: 10.1007/s40708-014-0002-y. [PMC free article] [PubMed] [Cross Ref]
17. Zong X, Lee J, Poplawsky AJ, Kim SG, Ye JC. Compressed sensing fMRI using gradient-recalled echo and EPI sequences. NeuroImage. 2014;92:312–321. doi: 10.1016/j.neuroimage.2014.01.045. [PMC free article] [PubMed] [Cross Ref]
18. Candes E, Romberg J, Tao T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Trans Inform Theory. 2006;52(2):489–509. doi: 10.1109/TIT.2005.862083. [Cross Ref]
19. Afonso M, Sanches JM. Image reconstruction under multiplicative speckle noise using total variation. Neurocomputing. 2015;150(Part A):200–213. doi: 10.1016/j.neucom.2014.08.073. [Cross Ref]
20. Shen Y, Li J, Zhu Z, Cao W, Song Y. Image reconstruction algorithm from compressed sensing measurements by dictionary learning. Neurocomputing. 2015;151(Part 3):1153–1162. doi: 10.1016/j.neucom.2014.06.082. [Cross Ref]
21. Eslahi N, Aghagolzadeh A, Andargoli SMH. Image/video compressive sensing recovery using joint adaptive sparsity measure. Neurocomputing. 2016;200:88–109. doi: 10.1016/j.neucom.2016.03.013. [Cross Ref]
22. Gao X, Jiang F, Liu S, Che W, Fan X, Zhao D. Hierarchical frame based spatial temporal recovery for video compressive sensing coding. Neurocomputing. 2016;174(Part A):404–412. doi: 10.1016/j.neucom.2015.07.110. [Cross Ref]
23. Candès EJ, Li X, Ma Y, Wright J. Robust principal component analysis? J ACM. 2011;58(3):11:1–11:37. doi: 10.1145/1970392.1970395. [Cross Ref]
24. Chandrasekaran V, Sanghavi S, Parrilo PA, Willsky AS (2010) Rank-sparsity incoherence for matrix decomposition. Technical report
25. Moore B, Nadakuditi R, Fessler J (2014) Improved robust PCA using low-rank denoising with optimal singular value shrinkage. In: Statistical Signal Processing (SSP), 2014 IEEE Workshop on, pp 13–16. doi:10.1109/SSP.2014.6884563
26. Jolliffe I. Principal component analysis. Hoboken: Wiley; 2005.
27. Scharf LL. The SVD and reduced rank signal processing. Signal Process. 1991;25(2):113–133. doi: 10.1016/0165-1684(91)90058-Q. [Cross Ref]
28. Tufts D, Shah A. Estimation of a signal waveform from noisy data using low-rank approximation to a data matrix. IEEE Trans Signal Process. 1993;41(4):1716–1721. doi: 10.1109/78.212753. [Cross Ref]
29. Drineas P, Kannan R, Mahoney MW. Fast monte carlo algorithms for matrices II: computing a low-rank approximation to a matrix. SIAM J Comput. 2006;36(1):158–183. doi: 10.1137/S0097539704442696. [Cross Ref]
30. Klema V, Laub A. The singular value decomposition: its computation and some applications. IEEE Trans Autom Control. 1980;25(2):164–176. doi: 10.1109/TAC.1980.1102314. [Cross Ref]
31. Nadakuditi R. Optshrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage. IEEE Trans Inf Theory. 2014;60(5):3002–3018. doi: 10.1109/TIT.2014.2311661. [Cross Ref]
32. Liang ZP (2007) Spatiotemporal imaging with partially separable functions. In: Biomedical Imaging: From Nano to Macro, 2007. ISBI 2007. 4th IEEE International Symposium on, pp 988–991. doi:10.1109/ISBI.2007.357020
33. Amaldi E, Kann V. On the approximability of minimizing nonzero variables or unsatisfied relations in linear systems. Theoret Comput Sci. 1998;209(12):237–260. doi: 10.1016/S0304-3975(97)00115-1. [Cross Ref]
34. Donoho D. Compressed sensing. IEEE Trans Inform Theory. 2006;52(4):1289–1306. doi: 10.1109/TIT.2006.871582. [Cross Ref]
35. Candes EJ, Recht B. Exact matrix completion via convex optimization. Found Comput Math. 2009;9(6):717–772. doi: 10.1007/s10208-009-9045-5. [Cross Ref]
36. Cai JF, Candès EJ, Shen Z. A singular value thresholding algorithm for matrix completion. SIAM J Optim. 2010;20(4):1956–1982. doi: 10.1137/080738970. [Cross Ref]
37. Gu S, Zhang L, Zuo W, Feng X (2014) Weighted nuclear norm minimization with application to image denoising. In: Computer Vision and Pattern Recognition (CVPR), 2014 IEEE Conference on, pp 2862–2869. doi:10.1109/CVPR.2014.366
38. Moran JM, Jolly E, Mitchell JP. Social-cognitive deficits in normal aging. J Neurosci. 2012;32(16):5553–5561. doi: 10.1523/JNEUROSCI.5511-11.2012. [PMC free article] [PubMed] [Cross Ref]
39. Zhang S, Block KT, Frahm J. Magnetic resonance imaging in real time: advances using radial flash. J Magn Reson Imaging. 2010;31(1):101–109. doi: 10.1002/jmri.21987. [PubMed] [Cross Ref]
40. Geethanath S, Reddy R, Konar AS, Imam S, Sundaresan R, Venkatesan R. Compressed sensing MRI: a review. Crit Rev Biomed Eng. 2013;41(3):183–204. doi: 10.1615/CritRevBiomedEng.2014008058. [PubMed] [Cross Ref]
41. Lingala SG, Hu Y, DiBella E, Jacob M. Accelerated dynamic MRI exploiting sparsity and low-rank structure: kt SLR. IEEE Trans Med Imaging. 2011;30(5):1042–1054. doi: 10.1109/TMI.2010.2100850. [PMC free article] [PubMed] [Cross Ref]
42. Graedel NN, McNab JA, Chiew M, Miller KL. Motion correction for functional MRI with three-dimensional hybrid radial-cartesian EPI. Magn Reson Med. 2016 [PubMed]
43. Hansen PC. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev. 1992;34(4):561–580. doi: 10.1137/1034115. [Cross Ref]
44. Lin T, Herrmann F. Compressed extrapolation. Geophysics. 2007;72:77–93. doi: 10.1190/1.2750716. [Cross Ref]
45. Lindquist MA. The statistical analysis of fMRI data. Stat Sci. 2008;23(4):439–464. doi: 10.1214/09-STS282. [Cross Ref]
46. Friston K, Ashburner J, Kiebel S, Nichols T, Penny WE. Statistical parametric mapping: the analysis of functional brain images. London: Academic Press; 2006.

Articles from Brain Informatics are provided here courtesy of Springer