PEAR: PEriodic And fixed Rank separation for fast fMRI

Purpose In functional MRI (fMRI), faster acquisition via undersampling of data can improve the spatial‐temporal resolution trade‐off and increase statistical robustness through increased degrees‐of‐freedom. High‐quality reconstruction of fMRI data from undersampled measurements requires proper modeling of the data. We present an fMRI reconstruction approach based on modeling the fMRI signal as a sum of periodic and fixed rank components, for improved reconstruction from undersampled measurements. Methods The proposed approach decomposes the fMRI signal into a component which has a fixed rank and a component consisting of a sum of periodic signals which is sparse in the temporal Fourier domain. Data reconstruction is performed by solving a constrained problem that enforces a fixed, moderate rank on one of the components, and a limited number of temporal frequencies on the other. Our approach is coined PEAR ‐ PEriodic And fixed Rank separation for fast fMRI. Results Experimental results include purely synthetic simulation, a simulation with real timecourses and retrospective undersampling of a real fMRI dataset. Evaluation was performed both quantitatively and visually versus ground truth, comparing PEAR to two additional recent methods for fMRI reconstruction from undersampled measurements. Results demonstrate PEAR's improvement in estimating the timecourses and activation maps versus the methods compared against at acceleration ratios of R = 8,10.66 (for simulated data) and R = 6.66,10 (for real data). Conclusions This paper presents PEAR, an undersampled fMRI reconstruction approach based on decomposing the fMRI signal to periodic and fixed rank components. PEAR results in reconstruction with higher fidelity than when using a fixed‐rank based model or a conventional Low‐rank + Sparse algorithm. We have shown that splitting the functional information between the components leads to better modeling of fMRI, over state‐of‐the‐art methods.


I. INTRODUCTION
Accelerating acquisition in functional MRI (fMRI) has gained significant attention in neuroimaging.
Accelerated fMRI can provide data at higher frame rates or sampling bandwidths, leading to higher temporal degrees of freedom 1 .This enables the use of more powerful and sophisticated analysis techniques 2 .
Alternatively, accelerated fMRI may be used to increase spatial resolution without sacrificing temporal fidelity, enabling time-resolved studies of the functional organization of the brain at finer scales, like cortical layers 3 or columns 4 .In resting state fMRI, where the goal is to estimate brain connectivity networks of the subject, accelerated data acquisition can improve the estimation of resting state networks (RSNs) 5 .
Numerous methods for accelerating acquisition of MRI data by exploiting its intrinsic structure and redundancy have been published.In clinical dynamic MRI (e.g.Cardiac MRI and Dynamic Contrast Enhanced (DCE) MRI), many methods are based on undersampling in the k-t space [6][7][8] .Since the introduction of compressed sensing (CS) 9,10 , accelerated CS-based methods for clinical dynamic MRI have also been developed 11,12 .While some of these methods have been adapted to fMRI [13][14][15][16][17][18][19] , the different nature of the fMRI signal (when compared to clinical dynamic MRI), e.g. its low variance of signal of interest and its limited spatial compressibility, limits the adoption of those implementations.
Recently, we introduced an approach for reconstruction of fMRI from undersampled measurements that is based on modeling fMRI as fixed-rank, i.e., k-t FASTER (FMRI Accelerated in Space-time by means of Truncation of Effective Rank) 20 .It aligns with the common analysis approaches in fMRI, which estimate limited numbers of spatial and temporal components from the data.It has been demonstrated that k-t FASTER provides higher quality of activation and resting state maps when compared with other methods for fMRI reconstruction from undersampled data.This method has been extended recently to include multi-channel coil sensitivity information and more flexible radial-Cartesian sampling 21 , providing additional encoding information and more incoherent sampling of k-space, resulting in robust recovery of task-based fMRI data at acceleration factors of 10 times or higher.
Several other approaches for reconstruction of fMRI from undersampled measurements have been recently suggested [22][23][24][25][26] .Aggarwal et al. examined enforcing a low-rank model and signal sparsity 22 .Others explored exploiting low rank and sparsity after the separation of fMRI into two components.This approach, known as Low-rank plus Sparse (or L+S) 27 , consists of modeling the data as a sum of two components, where low rank is enforced on one of them, and sparsity in some transform domain is enforced on the other.L+S has been examined for both clinical dynamic MRI 28 and fMRI [23][24][25]29 .
The common implementation of L+S for clinical dynamic MRI and fMRI consists of modeling the low rank component as a "background" image while the sparse component contains the dynamic information.
This approach leads to satisfactory results for clinical dynamic MRI applications where the dynamic signal is significantly above the noise level (e.g.MR angiography (MRA) and DCE-MRI) or periodic in the time domain (e.g.cardiac MRI).However, based on our experiments, its performance for fMRI, where the signal is often near the noise level and filtered by the hemodynamic response, is sub-optimal.
In this study, we examine a different separation of the data, where both components contain functional information.While most previous methods that combine low rank and sparsity are based on solving an unconstrained minimization problem by singular value soft-thresholding (SVT) 30 , in our approach we solve a constrained minimization problem based on truncating the singular values (a.k.a Truncated SVD or TSVD) 31 .Our approach forces one of the components to have a moderate fixed rank, and the other to be sparse in the temporal Fourier domain, leading to improved results compared to an SVT-based approach.
Reconstruction is performed via alternating minimization, that enforces the fixed-rank requirement and sparsity iteratively.
We call our approach PEAR: PEriodic And fixed Rank separation for fast fMRI.We examine reconstructions from undersampled data acquired using golden-angle radial sampling 32 , and correspondence to both time-course information (using General Linear modeling (GLM)) and spatial information (resting state network maps estimated via dual regression 33 ).Our experiments consist of a purely synthetic simulation, to show the concept of separation between the components, a synthetic simulation using real timecourses to examine correspondence of results to real and known time-courses, and retrospective sampling of a real resting state fMRI dataset, to examine resting state network recovery.We compare PEAR to k-t FASTER that uses a fixed-rank model only, and to a conventional, SVT-based L+S implementation.We also explore the contributions of the different components in PEAR.Based on our experiments, PEAR exhibits better estimation of the timecourses and resting state networks from undersampled data, compared to k-t FASTER and to L+S, using only 6.25% of the data in the simulations and 10% of the data in the retrospective sampling experiment.
The paper is organized as follows.Section II presents the proposed method for faster fMRI via separation of signals into periodic and fixed rank components.Section III describes experimental results.Section IV discusses theoretical aspects and implementation details of our method and Section V concludes by highlighting the key results.

II. METHOD
In MRI, data is acquired in the spatial Fourier domain (k-space).In dynamic MRI applications, such as cardiac MRI, MRA and DCE MRI, as well as in fMRI, the k-space of each temporal frame is acquired.
By undersampling k-space (i.e.taking only partial k-space measurements for each temporal frame), one can obtain higher frame rate, or alternatively, cover a greater extent in k-space, thereby increasing spatial resolution without decreasing temporal resolution.
In the problem of fMRI reconstruction from undersampled k-space, our goal is to recover the time series of acquired images.For simplicity, the time series is represented as a space-time matrix, X ∈ R N ×T where each column is a 3D temporal frame concatenated as a vector, N denotes the number of pixels in a single frame, and T denotes the number of frames in the time series.The measurement model, which takes into account that in most cases data is acquired using multiple coils is: where y is a vector of undersampled measurements and E is a general linear operator that maps a matrix to a vector.For acquisition with multiple receiver coils, E consists of multiplication by coil sensitivities followed by an undersampled Fourier transform.The vector z represents the measurement noise, modeled as complex Gaussian with zero mean.
Since y is generated via undersampling, proper reconstruction of X from y requires assumptions on X.
Relying on framework of CS, many methods that are based on sparsity of X in some transform domain were examined for dynamic MRI in general and for fMRI in particular.In our recent work, we considered modeling X as a fixed rank matrix, which aligns with the theory that X is composed of a relatively small number of spatially coherent temporal processes.The fixed-rank based approach for fMRI (i.e., k-t FASTER) solves the following minimization problem 20 : min where r is a fixed, moderate rank that ranges between 20 and 50 in our fMRI approach (but may be much lower in other MRI modalities).Unlike some other types of dynamic MRI that exhibit high variance of signals of interest, in fMRI valuable information may also be embedded in low variance components, slightly above the noise level.Consequently, method to retrieve information from higher dimensions is expected to provide better results for fMRI.
An approach that was applied initially for clinical dynamic MRI 28 , and has been examined recently for fMRI 23,24,29 , consists of modeling the dynamic sequence as a sum of two components.A low-rank component that represents mainly the background (denoted as the L component), and a component that contains the valuable signal.The latter is modeled as sparse in some transform domain, and denoted as the S component.This approach, known as L+S 27 , is based on solving the following unconstrained problem: min where L and S denote the low-rank and sparse components, • * and • 1 are the nuclear norm and the ℓ 1 norm, λ 1,2 are tuning parameters that control the weight given to each term in the optimization problem and the reconstructed space-time matrix is X = L + S. The linear transformation Ψ is a sparsifying transformation applied on S, depending on the specific dynamic MRI application.For MRA, Ψ may be chosen as an identity transform, whereas for cardiac MRI, which consists of periodic temporal structure, Ψ may be a temporal Fourier transform (i.e., applying a Fourier transform row-wise, independently on each of the rows of S) 28 .For fMRI, both types of transformations were examined: Otazo et el. 29considered the identity transform (although their decomposition is used for analysis rather than acceleration) and Singh et al. 23 examined the temporal Fourier transform.We note that L+S solves an unconstrained problem that does not explicitly enforce a fixed rank, and the solution is often based on SVT.By viewing the results of current implementations of L+S for fMRI 23,29 we found that in practice the resulting L component tends to have a very low rank and typically contains only background information, while the important functional information is in the S component.
Some fMRI analysis models suggest that while neural signals have strong band limited components, they also exist across the frequency spectrum 34 .Therefore, we consider including sparsity in the temporal spectrum, to capture the bandlimited assumption, in addition to a fixed rank representation that would be more suitable for broader-band signals.
In particular, we propose modeling the fMRI signal as a sum of a fixed rank component, which contains the high variance information, and a periodic component that captures the periodicity that is not captured in the fixed rank component.Thus, we model the fMRI data X as X = A + P, where A and P are the fixed rank and periodic components, respectively.We enforce a limited number of temporal periodic signals for P, by demanding sparsity in the temporal Fourier domain, and a fixed rank for A which contains the high variance signal.
To understand the rationale behind this modeling, Fig. 1   Fig. 1 Top: Timecourses (left) and amplitude spectrum (right) after removing DC of four selected pixels from a preprocessed resting-state fMRI dataset.Bottom left: spatial locations of selected pixels.Pixels #1-#3 showed high correspondence (|Z| > 6) with Default Mode Network and pixel #4 showed high correspondence with a visual network.The amplitude spectrum of the pixels shows that while pixels #1, #2 and #4 exhibit a limited number of peaks in the spectral domain (and therefore may be suitable for sparse modelling in the temporal Fourier domain) pixel #3 involves a wide range of spectral components.As a result, separation into a fixed-rank and periodic components as proposed by our method would allow better representation of those signals, compared to using a fixed rank component (k-t FASTER) or a periodic component only (L+S).
To obtain the separation into A and P components, we propose the following minimization problem: where C is the set of matrices with a fixed rank r (ranges between 20 and 50 in our fMRI approach) and F t is the temporal Fourier transform that applies a Fourier transform row-wise on each of the rows of P (where X = Â + P).
We solve (4) using alternating minimization (AM) 35 .In this approach, in each iteration we perform minimization with respect to one variable while keeping the other one fixed, and then switch between the variables.In our case, we start with an arbitrary initial point P 0 .For n ≥ 1 we iteratively compute: (5) We solve each sub-problem (5,6) via gradient projection 36 where the proximal gradient is used for the nondifferentiable ℓ 1 function in (6).A solution for ( 5) has been proposed by Goldfarb et al. 37 , a.k.aIterative Hard Thresholding with Matrix Shrinkage (IHT-MS).It consists of a gradient step for data consistency followed by a projection step onto the subspace C. The general step is: where R r (Q) is the projection onto the subspace C, defined as: R r the singular values of Q, and u i and v i are the singular vectors associated with σ i .The operator S µ (Q) is the singular value soft-thresholding operator, defined as: .., σ m ) and U and V are the left and right singular vectors associated with Σ.The parameter α is a step size, which controls the rate of convergence.The rationale behind applying S µ (•) before applying R r (•) can be explained by modelling the lower singular values, {σ i } m i=r+1 as representing noise, while {σ i } r i=1 represent data contaminated with additive noise.The subtraction of µ from {σ i } r i=1 is explained as removing noise from the higher singular values.Based on our experiments 20 , the selection of µ = c • σ r+1 where c is a fixed parameter, leads to good results.
To solve (6) we use the Iterative Shrinkage-Thresholding Algorithm (ISTA) 38,39 , whose details are given in Appendix A. The general step for the solution of ( 6) is: where Λ λ indicates the soft-thresholding operator with parameter λ, applied element-wise.Finally, by defining 7) and ( 8) become: 9) The convergence of an iterative solution based on (7) has been proven by Goldfarb et al. and To summarize, the major differences between L+S and the PEAR approach for fMRI are outlined below: • Thresholding mechanism and a fixed rank solution: The solution of the L+S problem given in (3) using the same AM approach used for solving (4) results in an algorithm that is different from Algorithm 1 in the thresholding mechanism of the singular values.While in Algorithm 1 we perform singular value soft-tresholding (SVT) with the value c • σ r+1 followed by truncating the r + 1...m singular values (TSVD), for the solution of (3) we perform SVT with value λ 1 , with no rank constraint.The issue of an SVT-based solution versus a TSVD-based solution has been studied in the literature previously 42 , and it has been shown that TSVD performs better for fixed rank problems 42 .
Based on our experience, fMRI can be considered as a fixed rank problem, since the number activation networks is relatively small and in many cases kept fixed for analysis.This statement also aligns with our experimental results for fMRI hereinafter.
• Separation of functional information between components: In conventional implementation of L+S for fMRI 23,24,29 , the L component is modeled as very low rank, and therefore contains background information and no functional information.In PEAR, functional information is split between the components.Consequently, recovered functional information is not limited to periodic signals only, and results are improved compared to L+S, as will be shown in the next section.Indeed, the solution of ( 4) can be approximated by solving (3) with appropriate selection of λ 1,2 values.However, PEAR enables enforcing a fixed rank (based on a priori knowledge of typical fMRI dimensionalities, often between 20-50) for a variety of datasets, obviating the need to examine a range of λ 1,2 values for each separate dataset.In addition, our experiments show that solving (4) provides better results for fMRI, when compared to solving (3) also for the case where λ 1,2 were carefully chosen for optimality.

III. EXPERIMENTAL RESULTS
To demonstrate the performance of PEAR compared to a well defined ground truth and additional algorithms for the reconstruction of fMRI from undersmpled measurements, we performed 3 types of experiments.In all experiments PEAR is compared to k-t FASTER 20 and L+S 23 (where sparsity is enforced in the temporal Fourier domain, as shown in ( 3)).
The first experiment is a simulation based on synthetically generated mixtures of periodic and aperiodic signals, and aims to compare the results of PEAR to the aforementioned methods and to examine how PEAR separates the signals into periodic and fixed rank components.The second experiment is an extension of the first experiment using realistic time-courses instead of purely synthetic ones.In the third experiment we examine the performance of PEAR for a retrospectively undersampled realistic 3D fMRI data sequence.
In all experiments, data is undersampled retrospectively, through a radial sampling approach using the golden-angle 32,43 and the NUFFT 44 package was used for forward and adjoint spatial Fourier transforms.
The output of each experiment is provided as z-statistics maps that reflect the degree to which each timecourse (in the case of experiments 1 and 2) or spatial regressor (in case of experiment 3) is expressed with a unique time-course in the data.Output maps were null-corrected using a Gaussian and Gamma mixture model 5 .
The parameters c = 0.7, α = 1 (for k-t FASTER) and α = 0.5 (for PEAR and L+S) were selected experimentally.In all cases, all time-points were initialized to the mean image calculated from all projections.All algorithms were run 100 iterations or until the minimum update between consecutive iterations was below 10 −4 .For k-t FASTER and PEAR we examined the parameter r in the range between 1 and 50, and experimentally used r = 32 for k-t FASTER, r = 27 for PEAR in experiments 1 and 2, and r = 20 for PEAR in experiment 3.In experiments 1 and 2 the parameters were tuned for optimal performance for each method (where optimal performance is evaluated by examining the z-stat maps), and in experiment 3 parameters were tuned for optimality on a training set and results were obtained using the same parameters for an unseen fMRI sequence.

III.A. Experiment 1: Purely synthetic simulation
In this experiment, we simulated a phantom consisting of 5 Regions of Interest (ROIs).Each ROI is spatially formed as a single letter from the letters "FMRIB", and contains one of 5 purely synthetic timecourses, generated as follows.The letters "F" and "I" were purely periodic timecourses where "F" contains a single frequency and "I" contains a mixture of three frequencies, "R" was a purely aperdioc timecourse, and "M" and "B" were a superposition of periodic and aperiodic timecourses.The phantom was added to a realistic background fMRI dataset, to form a 2D fMRI sequence with known functional timecourses, of size 64 × 64, with 512 time points.The timecourses and their spatial locations in the simulated image are shown in Fig. 2 (left and top).
Undersampling was carried out in the k-t space.We examined two undersampling ratios, first by taking 8 radial projections at each timepoint (corresponding to acceleration ratio of R=8 relative to a fully-sampled, maximally efficient Cartesian acquisition).We then repeated the experiment using only 4 radial projections at each timepoint (corresponding to R=16).This simulates one slice of a hybrid radial-Cartesian trajectory, which rotates an EPI trajectory within 3D k-space 21 .An additive white Gaussian noise with zero mean was added to the samples in the k-space domain to obtain SNR of 25dB.For PEAR, λ was examined in the range of 0.45-3.4and was selected as λ = 0.91 experimentally.For L+S, λ 1 was examined in the range of 1.1-3.4 and λ 2 was examined in the range of 0.45-3.4.These parameters were selected as λ 1 = 1.6 and λ 2 = 0.91 experimentally (the values for λ, λ 1,2 are provided after normalization with respect to the standard deviation of the data).
To examine the correspondence of the reconstruction with the ground-truth, we performed regression against the original timecourses using General Linear Model (GLM) 45 .Figure 3 shows the F-test results as null-corrected z-statistics maps 5 , for the ground-truth data (fully sampled image without the addition of noise), L+S, k-t FASTER and PEAR, for both R=8 and R=16.All maps are thresholded at |Z| > 4.3 and shown with color scale mapped between 4.3 < |Z| < 15.
It can be seen that PEAR provides the most reliable result, being the only method that almost perfectly recovers both "M" and "B" with minimum false positive errors at R=16.In addition, we see that L+S is unable to recover the aperiodic timecourse "R", as opposed to both k-t FASTER and PEAR thanks to their fixed-rank component.
An interesting analysis is the contribution of each component in PEAR. Figure 4 shows the GLM results for the A and P components of PEAR separately (for R=16), where the z-statistics maps are thresholded at |Z| > 4.3 and shown with color scale mapped between 4.3 < |Z| < 15.Note that the sum of the null-corrected z-statistics maps of A and P is not equal to the z-statistics map of PEAR, due to the null-correction applied for each map, that depends with each map's noise level.However, the separation of PEAR into periodic and fixed rank components is clearly demonstrated.The A component highly corresponds with the letter "R" which is a purely aperiodic timecourse, and with the letters "M" and "B" that include an aperiodic part.The P component highly corresponds to the letters "F" and "I" which are purely periodic timecourses, and to the letters "M" and "B" that include an periodic part.As demonstrated, this separation allows better modelling and leads to better recovery compared to k-t FASTER and L+S.Another analysis is presented in Fig. 5, where example portions of the mean timecourses from the five letter ROIs are shown for the ground truth, L+S, k-t FASTER and PEAR reconstruction results, including the timecourses for the A and P components of PEAR separately.The timecourses are shown in arbitrary units, to allow proper examination of their structure.It can be seen that as expected, L+S is limited in its ability to track the rapid changes that appear in the letter "R".In addition, the P component of PEAR indeed contains the periodic part of the signal, and therefore exhibits high correspondence with letters that are fully periodic ("F" and "I").
These simulations clearly demonstrate the expected behaviour of our proposed approach under conditions where signals include pure periodicity; however, these are not a realistic depiction of fMRI data, even in task conditions, since these signals are rarely strongly periodic.The following experiment examines the performance of the various algorithms for realistic timecourses.

III.B. Experiment 2: Simulation with realistic timecourses
In this experiment, we generated 5 realistic timecourses.First, we used a regression of a realistic fMRI dataset taken from the Human Connectome Project (HCP) database 46 against 15 canonical resting It can be seen that L+S is unable to recover the letter "R", due to its purely aperiodicity.While PEAR exhibits better results for R=8 and R=16 when compared to the other methods, the difference between k-t FASTER and PEAR is emphasized for R=16, where PEAR provides the most reliable result, with almost perfect recovery of the letters "M" and "B", at minimal ratio of false positive errors.It can be seen that the A component captures the letter "R" which is purely aperiodic, and the letters "M" and "B" that include an aperiodic part.The P component captures the letters "F" and "I" which represent periodic timecourses, and the letters "M" and "B" that include an periodic part.

PEAR: A component PEAR: P component
state network maps (RSNs) 33 , derived from high-dimensional group-level Independent Component Analysis (ICA) of resting state fMRI datasets.The regression result provided 15 timecourses, each one corresponds to a single RSN regressor.We pulled 5 timecourses and used them instead of the purely synthetic timecourses used in experiment 1.These timecourses were used for the same simulation of the letters FMRIB, rather than retained in the original network spatial maps, due to the ease of visually evaluating the output parameter maps.).The timecourses are shown in arbitrary units, to allow proper examination of their structure.It can be seen that L+S is limited in its ability to track the rapid changes that appear in the letter "R".In addition, the separation of PEAR into A and P component is clearly seen, as the P component exhibits high correspondence with letters that are fully periodic ("F" and "I"), and A exhibits high correspondence with the aperiodic letter, "R".
We repeated the same setting of experiment 1 (including the same SNR, sampling ratios, and selected parameters for each algorithm) where the only difference is the use of realistic timecourses instead of simulated ones.The realistic timecourses, including the corresponding regressors used for generation of the timecourses, are shown in Fig. 2 (bottom and left).
Figure 6 shows the General Linear Model (GLM) F-test 45 results as null-corrected z-statistics maps that were computed against the realistic time courses of all letters, for the ground-truth data (fully sampled image without the addition of noise), L+S, k-t FASTER and PEAR, for R=8 and R=16.All maps are thresholded at |Z| > 4 and shown with color scale mapped between 4 < |Z| < 15.
In correspondence with experiment 1, for realistic timecourses we see that for R=8, both k-t FASTER and PEAR provide similar results that outperform L+S.For R=16, we see that PEAR provides cleaner results, as can be seen mainly when comparing the recovery of the letters "F" and "B".It can also be seen that all methods are unable to recover "R" and "I" due to the high undersampling ratio.It can also be seen that "R" and "I" are the letters with the lowest Z values in the fully-sampled, ground truth data, due to their low energy in this experiment.

III.C. Experiment 3: Retrospective undersampling of real fMRI dataset
In this experiment, we examined kt-FASTER, L+S and PEAR for the scenario of undersampled 3D fMRI data, taken from the HCP database 46 .Data was registered to an MNI standard space with dimensions 91 × 109 × 91 and included 512 timepoints.We examined two sampling ratios.First, we used 15% (R=6.66) of the data by taking only 14 radial blades at each timepoint for each axial slice.In addition, we examined the scenario of using only 10% of the data (R=10).To keep memory requirements under control, reconstructions were performed independently for each 91x109 2D axial slice with 512 timepoints.
After reconstruction, all 2D reconstructed slices were stacked together to form a 3D image with 512 time point for further analysis.
After data reconstruction, we evaluated correspondence of the data to 15 canonical Resting State Networks (RSN) derived from high-dimensional group-level ICA of resting fMRI datasets from the HCP database 46 .Evaluation was done using dual regression 47 as follows.First, we performed spatial regression of the reconstructed dataset against the canonical maps (regressors) to extract the timecourses corresponding to each map.Then, we performed temporal regression of the dataset against the time series.The output is a set of z-statistic maps (one for each regressor) that reflect the degree to which each spatial regressor is expressed with a unique time-course in the data.
We compared the z-statistics maps from PEAR to those computed from the fully sampled data (ground truth) and from reconstructions using L+S and k-t FASTER methods.In this experiment, the parameters for each algorithm were tuned for a training sequence, and the results are evaluated using the same parameters for an unseen data sequence.For PEAR, λ was examined in the range between 0.25 and 2, and was selected as λ = 1.75 experimentally.For L+S, λ 1 was examined in the range between 0.25 and 1.75, and λ 2 was examined in the range between 0.13 and 1.75.These parameters were selected as λ 1 = 1.25 and λ 2 = 0.25 experimentally (the values for λ, λ 1,2 are provided after normalization with respect to the standard deviation of the data).
Figure 7 shows the Default Mode Network (DMN) regressor used for dual regression overlaid on the MNI atlas, as well as the z-stat dual regression outputs of the ground truth, kt-FASTER, L+S and PEAR for R=6.66.All maps are thresholded at |Z| > 3.3 and with color scale mapped between 3.3 < |Z| < 8.
The results for R=10, are shown in Fig. 8.The green ellipses in the images show regions where PEAR's activation pattern is the most similar one to the ground truth.It can also be seen that while both k-t FASTER and PEAR provide reliable results for both R=6.66 and R=10, L+S is does not provide satisfactory results at the higher acceleration ratio of R=10.
To provide measure for comparison between the methods, we computed the receiver operating characteristic (ROC) curves.This curve shows the performance of each method when compared to the ground truth (in terms of true positive ratio (TPR) against false positive ratio (FPR)) as the discrimination threshold varies.As a reference, we used the ground truth DMN z-stat map shown in Fig. 7 (thresholded at |Z| > 3.3).For the generation of ROC we computed the TPR and FPR for each DMN z-stat map for each method, as the threshold Z value ranges between 0 and 10.A common scalar measure for the performance of the algorithm is the area under the curve (often referred to as AUC). Figure 9 shows the ROC curves for kt-FASTER, L+S and PEAR, including the AUC computed for each curve, for R=6.66 and R=10.It can be seen that PEAR provides the most convex shape with the highest AUC in both cases.
In addition, the degraded performance of L+S for R=10 can clearly be seen by examining its ROC curve for R=10.To check the validity of this result for different RSN maps, we computed the AUC for all the 15 canonical RSN maps used in our dual regression process.The summary of the results is shown in Table I, where the method that provides the highest AUC for each RSN in marked in bold (for R=6.66 and R=10 separately).We see that while there are cases in which k-t FASTER or L+S outperform PEAR for some of the maps, PEAR provides the best AUC for the majority of the maps for both R=6.66 and R=10.The results of the training dataset for R=6.66 are also shown in the table, for completeness.
Finally, we examined the separation of PEAR into into A and P components, in terms of both time courses and spatial z-stat maps (for R=6.66).For this purpose, we first performed dual regression analysis for the A component and for the P component separately, to generate a z-stat map for each.Those maps, thresholded at |Z| > 3.3 and with color scale mapped between 3.3 < |Z| < 8, are shown in Fig. 10 and demonstrate that both A and P components contain functional activity.Then, we arbitrary selected a single pixel that exhibited high correspondence with DMN for both A and P (z-stat value of |Z| > 4.5).
Figure 11 shows the timecourses and amplitude spectra of the A and P components, for the selected pixel, where the spatial location of the pixel is shown at the bottom of the Figure .The timecourses and amplitude spectra show that the A component contains a wide range of frequencies, with some strong peaks in the spectrum for frequencies that have strong total energy.The P component contains a limited number of temporal frequencies and captures the low-energy periodicity that is not captured in A. This separation shows that a selection of a fixed, moderate rank for the A component, in addition to a demand for a L+S k-t FASTER PEAR  7).It can be seen that L+S provides very noisy results at this undersampling ratio.
limited number of temporal frequencies for the P component leads to the desired separation, which results in improved z-stat results shown earlier.

IV.A. Relation to previous works
As described in the Introduction, a few methods that have been published recently also focus on separation of fMRI into two components 23,25,26,29 .The main differences between our work and these prior works lie in the methodology and experiments.First, we enforce both components to contain functional information explicitly, by solving an unconstrained minimization problem that enforces the A component to have a fixed moderate rank, using a TSVD based solution.This is in contrast to L+S-based methods that practically enforce all the functional information to be contained in a single component, required to  includes the examination of the functional information that is contained in each of the components, is expanded compared to previous papers that deal with separation of fMRI into components; in particular, some of them show task-based MRI where the design is periodic or basic RSN analysis.

IV.B. Error measures and reproducibility
It is important to note that in the case of fMRI, conventional measures between the reconstructed datasets (e.g.MSE or correlations) may be misleading, as lower MSE does not necessarily mean that low variance functional information is preserved.Therefore, evaluation of results in this work is based the z-stat analysis of activation maps, which is the common tool used today for resting state fMRI analysis.
In addition, while prospective undersampling is possible (and has been carried out in our previous works 20 ), its evaluation has to be performed against connectivity maps taken from the literature or groupaveraged RSN spatial maps, leading to uncertainty in ensuring that subject-specific detail is retained.
Therefore, we focus on retrospectve undersampling, which allows accurate comparison against a well defined ground truth.It can be seen that the A component contains a wide range of frequencies, with some strong peaks in the spectrum for frequencies that have strong total energy.The P component contains a limited number of temporal frequencies and captures the periodicity that is not captured in the fixed-rank component, A. This separation shows that a selection of a fixed, moderate rank for the A component, in addition to a demand for a limited number of temporal frequencies for the P component leads to the desired separation, which results in improved z-stat results shown earlier.

IV.C. Limitations
This work focuses on resting state fMRI, which is a branch in fMRI research that deals with mapping brain connectivity based on an fMRI experiment that does not involve stimulation.In contrast, taskbased fMRI involves known manipulations of brain activity (a.k.a task-based fMRI), although some of the properties of task-based fMRI data are very similar to resting fMRI (i.e. the low variance signal based on the BOLD effect).Since many other reconstruction techniques place strong assumptions on brain activity, e.g. that it is periodic, highly reproducible or smooth in time, PEAR does not place these assumptions and therefore is expected to provide reliable results also for task-base fMRI.However, the analysis of task-based fMRI with PEAR is reserved for future research.
In addition, the reconstruction time of PEAR for a single realistic 2D axial slice in the dimensions described in our experimental results (109 × 91 with 512 time points) is approximately 15 minutes using MATLAB running on a single machine with 3.2GHz CPU.To allow recovery of a multi-slice image, we used cluster-based computing that processed all 91 slices in parallel and provided a 3D volume approximately at the same running time.While the computation time and the need for cluster-based computing are drawbacks of all iterative methods used to reconstruct fMRI sequences examined in this paper (k-t FASTER and L+S), it is performed offline, while the subject is no longer in the scanner and does not require expensive MRI resources.In addition, we are examining approches to accelerate the reconstruction process, mainly by improving the NUFFT, using methods with lower computational complexity 48 and implementing algorithms in a real-time programming environment.

V. CONCLUSIONS
This paper presents PEAR, an under-sampled fMRI reconstruction approach based on separating the fMRI signal to periodic and fixed-rank components.The higher accelaration ratio offered by PEAR results in reconstruction with higher fidelity than when using a fixed-rank based model or a conventional L+S algorithm.We have shown that splitting the functional information between the A and P components, by solving a constrained problem that enforces a fixed, moderate rank for the A component, leads to better modeling for fMRI, due to the unique nature of the fMRI signal.Future work will focus on extending this work to task-based fMRI using both retrospective and prospective sampling.
shows 4 timecourses of arbitrarily selected pixels in a resting state fMRI dataset, that exhibited high correspondence (|Z| > 6) with regressors that represent a Default Mode Network (DMN) map or a visual network map after a dual-regression against those network maps.The signals in the figure were extracted from the fMRI sequence after conventional fMRI pre-processing (including skull stripping, motion correction and slice timing correction), and the spatial locations of the selected pixels are shown at the bottom left of the figure.The timecourses are shown in the time domain (top left) and in the temporal Fourier domain, after removing the DC component (top right).It can be seen that while three of the timecourses exhibit peaks in the spectral domain and are suitable for sparse representation in the temporal Fourier domain, one of the timecourses exhibits a broad spectrum in the temporal Fourier domain.As a result, separation into fixed rank and periodic components allows better representation of signals (compared to fixed-rank only or periodic component only) as it is expected to capture both broad-band and band-limited temporal spectra.

Fig. 2
Fig.2Left: Spatial locations of the ROIs used in experiments 1 and 2, each ROI is formed as single letter and contains a single timecourse.Top: The timecourses used for each ROI in experiment 1. "F" and "I" are purely periodic timecourses where "F" contains a single frequency and "I" contains a mixture of three frequencies, "R" is a purely aperdioc timecourse, and "M" and "B" are superposition of periodic and aperiodic timecourses.Bottom: Timecourses used in experiment 2 for each simulated ROI (letter).

Fig. 3
Fig. 3 Experiment 1: GLM F-test results of L+S, k-t FASTER and PEAR for purely synthetic simulation, for R=8 (top) and R=16 (bottom).The z-stat map of the ground truth is also shown (left).All maps are thresholded at |Z| > 4.3 and with color scale mapped between 4.3 < |Z| < 15.It can be seen that L+S is unable to recover the letter "R", due to its purely aperiodicity.While PEAR exhibits better results for R=8 and R=16 when compared to the other methods, the difference between k-t FASTER and PEAR is emphasized for R=16, where PEAR provides the most reliable result, with almost perfect recovery of the letters "M" and "B", at minimal ratio of false positive errors.

Fig. 4
Fig.4Experiment 1: GLM F-test results of A and P components of PEAR for R=16.Maps are thresholded at |Z| > 4.3 and with color scale mapped between 4.3 < |Z| < 15.It can be seen that the A component captures the letter "R" which is purely aperiodic, and the letters "M" and "B" that include an aperiodic part.The P component captures the letters "F" and "I" which represent periodic timecourses, and the letters "M" and "B" that include an periodic part.

Fig. 5
Fig.5Example portions of the mean timecourses of the five letter ROIs are shown for the ground truth, L+S, k-t FASTER and PEAR reconstruction results, including the timecourses for the A and P components of PEAR separately (Experiment 1, R=16).The timecourses are shown in arbitrary units, to allow proper examination of their structure.It can be seen that L+S is limited in its ability to track the rapid changes that appear in the letter "R".In addition, the separation of PEAR into A and P component is clearly seen, as the P component exhibits high correspondence with letters that are fully periodic ("F" and "I"), and A exhibits high correspondence with the aperiodic letter, "R".

Fig. 6
Fig. 6 Experiment 2: GLM F-test results of ground truth, L+S, k-t FASTER and PEAR for simulation with realistic timecoureses, for R=8 (top)and R=16 (bottom).The z-stat map of the ground truth is also shown (left).All maps are thresholded at |Z| > 4 and with color scale mapped between 4 < |Z| < 15.It can be seen that although "R" and "I" are almost irrecoverable and PEAR and k-t FASTER provide similar results for R=8, PEAR provides better results for R=16 with minimal ratio of false positive errors.

Fig. 7
Fig. 7 Experiment 3: Retrospective sampling of realistic fMRI dataset.Left: the regressor used for dual regression overlaid on the MNI template, and the dual regression results of the ground truth.Right: Dual regression results of L+S, k-t FASTER and PEAR obtained at undersampling ratio of 15% (R=6.66).All maps are thresholded at |Z| > 3.3 and with color scale mapped between 3.3 < |Z| < 8.The green ellipses illustrate regions where the activation pattern detected by PEAR is the most similar one to the ground truth.

Fig. 8
Fig. 8 Experiment 3: Retrospective sampling of realistic fMRI dataset.Dual regression results of L+S, k-t FASTER and PEAR obtained at undersampling ratio of 10% (R=10) .All maps are thresholded at |Z| > 3.3 and with color scale mapped between 3.3 < |Z| < 8.The green ellipses illustrate regions where the activation pattern detected by PEAR is the most similar one to the ground truth (shown in Fig.7).It can be seen that L+S provides very noisy results at this undersampling ratio.

Fig. 10 Z
Fig. 10 Z-stat maps of A and P components of PEAR for R=6.66.Maps are thresholded at |Z| > 3.3 and with color scale mapped between 3.3 < |Z| < 8.It can clearly be seen that both components contain functional information and present activation regions in locations that correspond to similar activation regions in the ground truth.

Fig.
Fig.11Top: Timecourses and amplitude spectra of a pixel resulted with |Z| > 4.5 for DMN z-stat map.Timecourses and amplitude spectra of ground truth (top row) and A and P components of PEAR result (for R=6.66) (second row) are shown.Values are shown in arbitrary units for better view.The spatial location of the selected pixel on an axial slice is also shown (left).It can be seen that the A component contains a wide range of frequencies, with some strong peaks in the spectrum for frequencies that have strong total energy.The P component contains a limited number of temporal frequencies and captures the periodicity that is not captured in the fixed-rank component, A. This separation shows that a selection of a fixed, moderate rank for the A component, in addition to a demand for a limited number of temporal frequencies for the P component leads to the desired separation, which results in improved z-stat results shown earlier.
the convergence of an iterative solution based on (8) is well studied in the literature38,40.Therefore, based on Csiszar and Tusndy 41 , the convergence of our proposed AM framework is guaranteed.The proposed algorithm is coined PEriodic And fixed Rank separation for fast fMRI (PEAR) and is summarized in Algorithm 1, where SVD represents the singular value decomposition and {σ i } r+1 i=1 are the singular values of the matrix X n−1 − P n−1 in descending order.The operator F H t is the conjugate temporal Fourier transform, and E H is the conjugate transpose of E.
t Predefined rank: r, Soft-shrinkage parameter: c Tuning constant: λ, Step size: α Iteration limit: N Output: Estimated fMRI time-series: X Initialize: P Fig. 9 Receiver operating characteristic (ROC) curve for experiment 3. Performance of L+S, kt-FASTER and PEAR are shown.The numbers in brackets indicate the area under the curve (AUC) for each method.The ROC results are generated for the ground truth map thresholded at |Z| > 3.3 serving as a reference.It can be seen that PEAR provides the highest AUC for DMN, for both R=6.66 and R=10.Table I Area under ROC for 15 RSN maps.The bold values indicate the method with the highest value in each line, for R=6.66 and R=10 separately.It can be seen that PEAR provides the best performance for most of the maps, for both R=6.66 and R=10.Training dataset results for R=6.66 are also shown for completeness.
be sparse in some transform domain.As a result, the L+S solution might be sub-optimal in reconstruction of signals that are neither periodic, nor strong enough to be captured in the low-rank component.Second, the experimental part of this paper presents a thorough validation of the suggested approach (compared to other existing methods), based on realistic nonuniform undersamping, and examines the spatial activation of resting state networks with broad spectrum characteristics.This analysis, which also A component P component