dataesr/scientific-paragraphs-categorization huggingface dataset:
g. , when more than 50 ICA components are estimated. In both cases, if multiple subjects are to be co-analysed (e. g. , in order to compare the restingstate networks between subjects), most researchers begin by carrying out a (low-or high-dimensional) ICA on the group dataset as a whole. The groupwise ICA components can then be mapped back onto individual subjects in order to allow for cross-subject network comparisons. For such data-driven decompositions of multi-subject datasets, the first (and most computationally intensive) stage in the analysis is normally to reduce the entire dataset down to a set of "group-average" spatial eigenvectors, using principal component analysis (PCA, based on singular value decomposition or SVD). The most natural approach here is to temporally demean and concatenate all subjects' datasets, and apply PCAeffectively treating all the data as if it were one single huge dataset.
For an n-dimensional group-ICA, the resulting n strongest spatial eigenvectors can then be fed into the core ICA unmixing algorithm, which identifies the set of group-average spatiallyindependent components. This last stage is not in general computationally burdensome, because the data now comprises at most a few hundred spatial eigenvectorsthe dataset has been reduced from the size of (voxels × timepoints × subjects) to (voxels × n). However, as increasingly large datasets/databases are being generated, with large numbers of timepoints and voxels in each subject's dataset, and very large numbers of subjects being combined, there are major computational challenges associated with group-level datadriven analyses. Current approaches for combining all subjects' datasets, and applying basic dimensionality reduction, cannot be run using the computational facilities available to most researchers, primarily because of the very large memory requirements.
There are two mathematically equivalent ways to estimate the PCA from the temporally-concatenated data from all subjects: a) estimate for each subject the voxels × voxels covariance matrix of temporal correlations (this is very large, as the number of voxels is large), and then average over subjects; or, b) estimate the (subjects × timepoints) × (subjects × timepoints) covariance matrix of spatial correlations, which will be very large if the number of subjects or number of individual-session timepoints is large. Existing approaches to this computational problem include the initial reduction of each subject's dataset to a (normally comparatively small) number of spatial eigenvectors, before these are then combined across subjects, and a final group-wise PCA computed  (Calhoun et al. , 2001).
Both methods' memory requirements do not increase with increasing numbers of subjects analysed, and overall execution time scales linearly with the number of subjects. We validate the accuracy and computational effectiveness of the approaches using rich simulations, comparing against alternative approaches, with favourable results. We also validate against an extremely large rfMRI dataset from the Human Connectome Project (HCP), utilising a 24-core compute server with 1. 25 TB of RAM to calculate the "gold standard" empirical result, achieving an accuracy of greater than 99. 99%. The first of the new approaches for large group-level PCA (referred to as MIGP -MELODIC's Incremental Group-PCA), was developed specifically with rfMRI in mind, and the validations described below use data that attempts to match the characteristics (dimensions, intrinsic rank, etc.
) of rfMRI data. Since the original development and validation of MIGP, we have become aware of closely related approaches being developed in computer science, where thorough theoretical investigations of the mathematical properties have been reported (with very positive conclusions regarding the accuracy of these approaches)  (Baker et al. , 2012; Řehůřek, 2010). The HCP has already started disseminating group-average "dense connectomes" (full voxelwise/vertexwise correlation matrices) calculated using MIGP, which is now released as part of the MELODIC ICA tool in FSL (FMRIB Software Library, www. fmrib. ox. ac. uk/fsl). The second, distinct, approach (referred to as SMIG -Small-Memory Iterated Group-PCA) builds on top of the original method proposed in  Hyvärinen and Smith (2012) , achieving increased accuracy compared with the original approach, by iterating the main estimation of group-average eigenvectors several times.
While both MIGP and SMIG have very low-memory requirements, and are accurate, MIGP may be more time-efficient if file I/O dominates total compute time, and SMIG would be more time-efficient otherwise. These approaches will hopefully also be of value for analysing large datasets from other modalities (and branches of science), particularly in this era of "Big Data"; however, here we have specifically tailored our detailed simulations to the characteristics of resting-state fMRI data. Note that we have specified resting-state specifically rather than also task-fMRI, because data reduction in the latter case is more straightforward computationally for large studiesin the cases where one does not have to account for the random temporal phase across datasets.
Method Background Each subject's dataset comprising t timepoints and v voxels can be represented as a 2-dimensional space-time matrix Y (t × v). We assume that the data has already been preprocessed to remove artefacts and align the data spatially into a standard space (co-ordinate system), so that the voxels are anatomically compatible across all subjects. We also assume that each voxel's timeseries has been demeaned (and potentially variance-normalised  (Beckmann and Smith, 2004) ). If we were to carry out a single-subject n-dimensional spatial-ICA, we would first apply a PCA, or equivalently the singular value decomposition, for example, representing the data as: Y tÂv ð Þ ≈U tÂn ð Þ Â S nÂn ð Þ Â V 0 nÂv ð Þ ; ð1Þ where n is typically much smaller than t, U is the set of temporal eigenvectors, V is the set of spatial eigenvectors, and the eigenvalues (components' strengths) lie on the diagonal of S (the above assumes that the t-n weakest components have already been discarded).
ICA is then applied to the matrix V, estimating a new set of spatial maps, which are linear combinations of the maps in V, and which are maximally independent of each other. In the case of a multi-subject dataset, a natural way to generate a group-averaged set of spatial eigenvectors (to feed into analyses such as group-average parcellation or group-ICA) is to temporally concatenate all s subjects' datasets, and apply PCA + ICA as described above, or, equivalently, to average the v × v matrices of temporal covariances across all subjects. The resulting PCA-based approximation will then be the same as above, but now with n × s "timepoints" in the temporally concatenated data. Unfortunately, with large datasets and large numbers of subjects, it becomes impractical to form this full concatenated dataset and then run a full PCA; memory limitations and/or computational time become prohibitive.
Various approaches have been suggested previously for large multi-subject fMRI datasets, and we now describe the most relevant prior work. Calhoun et al. (2001)  suggested approximating temporal concatenation by first reducing each subject's dataset to the top m spatial eigenvectors, then concatenating these eigenvectors (estimated separately for each subject) across all subjects, before running a final PCA to further reduce this dataset to the top n eigenvectors, to be fed into ICA. Typically m = n, meaning that each subject is reduced to the same number of PCA components as the final group-average ICA will estimate. Although using a small value for m limits the memory requirements for these operations, the data size does scale linearly with the number of subjects, which can eventually become impractically large.
MIGP -MELODIC's Incremental Group-PCA MIGP is an incremental approach that aims to provide a very close approximation to full temporal concatenation followed by PCA, but without the large memory requirements. The high accuracy is achieved by never reducing individual subjects' datasets to small numbers of PCA components. The incremental approach keeps an "internal" PCA space of m weighted spatial eigenvectors, where m is typically larger than the number of timepoints in each individual dataset. By "weighted", we mean that the eigenvalues (component strengths) are incorporated into the matrix of spatial eigenvectors. The final set of m components, representing the group-average (or temporally concatenated) PCA output, can then be reduced to the desired dimensionality n by simply keeping the top n components, and, if required, discarding the weightings (eigenvector amplitudes).
We now describe the approach in more detail. Start by temporally concatenating a number of subjects' datasets such that the number of combined timepoints is larger than m. Typically we might set m = 2 t, and begin by concatenating 2-3 datasets. This dataset is then fed into m-dimensional PCA, and the matrix W (m × v) = S (m × m) × V′ (m × v) is kept. Note that this is not just the set of (unweighted) eigenvectors; each is weighted by its respective eigenvalue, which is an important distinctionthe eigenvalues, characterizing the strength of an effect, are here used to scale the eigenmaps (spatial eigenvectors) so that the strength is retained in the reduced subspace. W becomes a running estimate of the final group-average set of spatial eigenvectors, and can be considered a pseudo-timeseries matrix of m "timepoints" and v voxels.
For each additional subject's dataset in turn, we now incrementally update W, by concatenating W with each dataset Y i , and applying PCA to reduce back to an updated W, again with just the top m components kept (see Appendix A). Because the eigenvalue weights are kept in W, this approach automatically achieves the "right" balance (relative weighting) between W and each new datasetwhere "right" means that the relative weighting between the running estimate and new datasets matches what would have happened in the fulltemporal-concatenation scenario. In other words, this retains the overall variance of each batch of data. 1  A computationally efficient approach for estimating the top m weighted eigenvectors is to first estimate the "timepoints" × "timepoints" covariance matrix, apply eigenvalue decomposition to extract the top m temporal eigenvectors (e.
e. , the pure mean of all raw timeseries matrices from all subjects). All subjects' rotated data matrices can then be averaged, and PCA applied, without ever needing to form large concatenated data matrices of timeseries or eigenvectorsalthough two "passes through" all of the original datasets is necessary. Building on this, and noting that the initial group-average data matrix will not have intrinsically "high resting-state correlation CNR" (due to the random phase of RSN timeseries in different datasets being averaged), we have developed an iterated approach to the estimation of the group-averaged spatial eigenvectors. We now describe the mathematical justification for this method, showing that the iterations converge towards the same output as PCA applied to full temporal concatenation.
To initialize the algorithm, we take the average over the Y i , W← 1 s X i Y i ð4Þ This implies choosing the first-stage dimension m to be equal to the number of timepoints (number of rows in Y i ). If we choose to use a smaller m, we can simply reduce the dimension of this initial value by ordinary PCA, without changing anything in what follows. 1 It is easy to show that the covariance of concatenated datasets equals the sum of the covariances of the individual datasets (even if they have different lengths and overall variance). Given the decomposition Y = USV′ = UW (where W is the weighted spatial eigenvectors), the (voxels × voxels) covariance is W′U′UW = W′W. Hence it is clear that when taking multiple (previously concatenated) datasets and concatenating these with a new single dataset, the concatenation followed by covariance estimation is equivalent to summing the two separate W′W estimates (one for the previously combined multiple datasets and one for the new dataset)hence no change in the relative weighting of the two separate covariances is necessary in order for the estimate of covariance to be equivalent to the temporal concatenation of a new subject onto a previously concatenated set of subjects.
Using group-PCA output to generate the "dense connectome" It is now trivial to show that one can easily estimate the large groupaverage "dense connectome" ("space" by "space" matrix of voxel-level temporal correlations) from the PCA output, should this be desiredstarting conceptually from the fully-temporally-concatenated dataset Y. Forming the dense connectome is not necessary if one only wishes to run group-level ICA, because normally all that is passed onto spatial-ICA, after the initial PCA, is the group-PCA spatial eigenvectors, as described above. However, if we do want to estimate the groupaverage v × v covariance matrix (dense connectome), we have: covariance Y ð Þ ¼ Y 0 Y≈VS 0 U 0 USV 0 ¼ VS 0 SV 0 ¼ W 0 W: ð9Þ where W can be the weighted spatial eigenvectors output by (e. g. ) MIGP.
Empirical evaluations We compare several group-PCA approaches: • TemporalConcatfull temporal concatenation of all subjects' datasets, followed by PCA. • MIGPas described above. • SMIGas described above. We need to specify m (the internal dimensionality retained) and a (the number of iterations). • GIFTthe method of  Calhoun et al. (2001) , concatenating withinsubject PCA outputs across subjects, followed by a final PCA. We also test GIFT with an internal dimensionality m double that of the final output dimensionality; this keeps more subject-specific detail than is normally done by default, at the cost of increased RAM requirement. • MeanProjectionan approach utilised in old MELODIC software versions, that projects individual datasets onto a PCA-reduced version of the group-average-data, concatenates the results across subjects, and then reduces down with a final group-level PCA.
RAM requirementsdifferent methods We first estimated the RAM required by several group-PCA methods, as a function of number of timepoints, voxels, subjects and estimated dimensionality. For each method: the size of the largest "timeseries" data matrix formed is estimated (in some cases, for example, this might be voxels × timepoints, in others, voxels × timepoints × subjects, and in others, voxels × subjects × dimensions). Additionally, the largest data covariance formed is estimated (either voxels × voxels, or timepoints × timepoints, depending on the method). Simple testing in MATLAB shows that running the most efficient SVD on a large data covariance matrix requires approximately double the total RAM than is needed to just store the covariance matrix. For our purposes then, we plot whichever is the larger of the two quantities; the largest timeseries data matrix formed, and the doubled covariance matrix size.
Fig. 1  shows these maximum-RAM estimates for a range of methods, under a number of different study scenarios. We use approximate/typical values for 4 study scenarios: typical small imaging studies; 1000 datasets from the "thousand functional connectomes" (KFC) data  (Biswal et al. , 2010) ; 1200 subjects (that will eventually be available) from the Human Connectome Project  (Van Essen et al. , 2013) , with a large number of timepoints per subject; and 100,000 subjects' datasets that will eventually be acquired by the UK Biobank Imaging study (www. ukbiobank. ac. uk), if it runs to completion as hoped. The voxel counts are either: 25,000 (number of brain voxels in MNI standard space, working at the rather crude resolution of 4 mm); 200,000 (2 mm MNI-space brain voxels); and 100,000 (the approximate number of grayordinates in the HCP standard co-ordinate system, with approximate 2 mm spacing between surface vertices and subcortical voxels).
Full temporal concatenation does not continue to rise greatly for the larger datasets, as might have been expected; this is because full temporal concatenation is mathematically equivalent to an approach based around averaging the (very large) within-subject voxels × voxels covariance matrices across all subjects. Therefore, once the dataset becomes large enough, this latter approach becomes the more RAM-efficient option, and hence the resulting RAM requirements stop increasing with the number of subjects. It can be seen that for none of these scenarios do the MIGP/SMIG methods require large amounts of RAMindeed never more than around 8 GB.
Simulationsmethods Our primary evaluations were based on simulated datasets containing a hierarchy of effects, including structured signal and artefact in the dataset for each simulated subject, and subject variability (including outlier effects). We generated several distinct simulations, spanning a range of different scenarios as described below. The core of the simulation is to randomly generate a number of (group-level ground truth) spatial maps (that will be modulated by randomly distinct timecourses), which in the simplest simulations will be exactly the same sets of maps for all subjects. The maps are randomly generated and somewhat sparse, with a unity-standard-deviation Gaussian random distribution embedded on top of values that are 0 in most voxels, and 5 in a minority of voxels.
e. , we have a single-subject outlier), to having multiple subjects in the second group. Individual subjects' datasets are then created on the basis of the "group-level" maps defined above. Each subject within a given group has spatial maps with controllable similarity to the group-maps defined above; each map for each subject has a controllable amount of subjectspecific variation added, and is then modulated by a single random timecourse (taking the outer product of the spatial map and a random  timecourse). All components are then added to form the full t × v data matrix. The subject-specific randomness can be added both to the spatial map (both its shape and amplitude) before temporal modulation, and after (the former therefore representing subject-specific alterations to the spatial map, compared with the group ground-truth, and the latter representing unstructured measurement white-noise).
Simulationsresults We now present results from a range of dataset scenarios, varying dataset dimensions, as well as the relative sizes of various forms of noise and subject variability. The full simulation parameter definitions are listed within each figure. In addition to the range of SMIG parameters reported on below, we also tested the following parameter settings; however, because they resulted in no (or very little) improved performance compared with the settings reported on, we do not include any further reporting on these. 1) With an internal (working) dimensionality m smaller than the number of original session timepoints t, we found virtually no accuracy advantage in raising the number of iterations from 10 to 25, and so do not report on more than 10 iterations. 2) With m = t, there was never a significant improvement when moving from 5 to 10 iterations, so we do not report on higher than 5.
Variations in subject variability and white noise (no artefacts) Fig. 2  shows results from relatively "high CNR, low variability" datasets (30 subjects, 200 timepoints). The first test has relatively low white noise and a low subject-variability of 0. 1 (meaning the extent to which the underlying common spatial maps vary across subjects). All methods perform well, being close to 100% accuracy in terms of both TPR (finding the correct common maps) and 1-FPR (only estimating the common maps and nothing extra). The second test increases the subject-variability to 0. 3, which might be considered to be relatively high (although high-resolution fMRI data with no spatial smoothing in the preprocessing might be expected to exceed this). All methods are still doing well, with around ~94% accuracy.
Interestingly, GIFT and MeanProjection perform very slightly better than the other methods. In the case of GIFT, this might be because of the potential "2-stage denoising" effect that could happen when the dataset exactly matches the assumptions and parameters being used by the method; given that GIFT is internally keeping exactly the correct number of components (here 10) in each within-subject PCA, and because this simulation has no subject-specific artefacts, then one expects that the estimated set of components for each subject should be close to perfect, achieving within-subject denoising effectively. Then, at the group level, when combining across the sets of withinsubject components, the higher level of subject variability can be more effectively dealt with, given that within-subject noise has been suppressed.
Interactions between subject-specific artefacts and estimated dimensionality Fig. 3  shows the effects of adding in subject-specific artefact components, and how the estimated dimensionality interacts with this. There are 10 common ("good") components in each subject's dataset, as well as 30 subject-specific "artefact" components (a realistic ratio  (Griffanti et al. , 2014) ). The difference between the 3 tests is the number of final estimated components, varying from being less than the true number (first test), equal (second test), and greater (third test). Not surprisingly, the TPR improves significantly when estimated dimensionality is raised, and 1-FPR falls. The performance of the best of the methods is now lower than before, with these more challenging scenarios, with iterated SMIG and MIGP reaching ~63-88% accuracies (depending on the exact scenario).
Very large datasets Fig. 6  shows tests with much larger datasets, with increased numbers of subjects or timepoints, and varying true and estimated dimensionality. The first and second tests have 50 (respectively, 200) subjects and 200 timepoints, with true and estimated dimensionality of 70 (a dimensionality that would most likely be too high to estimate robustly for individual subjects, but hopefully estimable at the group-level). GIFT with m = 2n, MIGP and iterated SMIG perform well, with a very slight difference between their best performances; here, the 10-iterations of m = t SMIG is 1-2% more accurate than the other options, and equal to TemporalConcat. The reason that here m = t, rather than m = 3n as in previous tests, is that here 3n N t. Default GIFT and MeanProjection perform less well than the other methods.
Subject ordering effects in MIGP Although there is no bias in overall eigenvector weighting as new subjects are added into the MIGP calculations, there still might be a small "ordering effect", for example, if eigenvectors found to be present early on (in the incremental estimations) may boost patterns that randomly match these in later datasets. Hence in this evaluation we simulated two fairly different groups of subjects. In addition to the standard MIGP analysis (which randomises the order that subjects are processed, and hence should avoid bias between different subject subgroups caused by any order effect), we also ran MIGP twice without randomisation. In one case ("forwards"), we first fed all subjects from group one into MIGP, followed by group two; in the other MIGP analysis ("backwards") we fed in group two first.
We also applied these 3 approaches using 5 iterations of the full algorithm, with each iteration initialised from the output of the previous one, rather than by the first subject. Fig. 7  shows the results from this evaluation. This shows that without order randomisation, some order bias can occur when different groups of subjects are processed, in group-organised ordering; small but significant biases can be seen in the results for "forwards" and "backwards" versions of MIGP, most strongly with respect to the accuracy of estimating group 2 in the "forwards" case (where group 1 is processed first). With the default MIGP approach of subject order randomisation, and only a single pass through the entire dataset, accuracy is very close to full temporal concatenation (which processes all subjects simultaneously).