dataesr/scientific-paragraphs-categorization huggingface dataset:
First, by looking across many datasets, we have learnt that encoding subspaces can sometimes be highly non-orthogonal (Figures  3-6e ) and hence not demixable under orthogonality constraints. Second, by not imposing orthogonality, we can easier identify components that are truly orthogonal. Third, removing the orthogonality constraint allowed us to obtain a simple analytical solution in terms of singular value decompositions (see Materials and methods) and hence to avoid local minima, convergence issues, and any additional optimization-related hyperparameters. To demonstrate these advantages, we ran the algorithm of  Brendel et al. (2011) , dPCA-2011, on all our datasets. The resulting components were similar to the components presented here, with the amount of variance captured by the first 15 components being very close; but the achieved demixing was worse.
For each component we defined a demixing index (see Materials and methods) that is equal to 1 if the component is perfectly demixed. For all datasets, these indices were significantly Table  1. Demixed PCA in comparison with existing methods. Columns: 'Signif. ' refers to the method of counting significantly tuned cells, as shown in Figure  1c-e. TDR refers to the 'targeted dimensionality reduction' of  Mante et al. (2013)  shown in Figure  1f-h. LDA stands for linear discriminant analysis, but this column applies to any classification method (e. g. support vector machine, ordinal logistic regression, etc. ). All classification methods can be used to summarize population tuning via a time-dependent classification accuracy (e. g. Meyers et al. , 2012). PCA stands for principal component analysis, as shown in Figure  1i-k.
, 2013). As mentioned above, the original application of this method yields only one component per task parameter and ignores the condition-independent components. While TDR can be extended in various ways to yield more components, no principled way of doing it has been suggested so far. Comparison of dPCA with TDR on our datasets shows that dPCA demixes the taskparameter dependencies better than TDR (see Figure  14  in the Materials and methods). For an in-depth discussion of the relationship between dPCA and LDA/MANOVA, we refer the reader to the Methods. Briefly, LDA is a one-way technique, meaning that only one parameter (class id) is associated with each data point. Therefore, LDA cannot directly be applied to the demixing problem. While LDA could be generalized to deal with several parameters in a systematic way, such a generalization has not been used for dimensionality reduction of neural data and does not have an established name in the statistical literature (we call it factorial LDA).
(2013) , in order to fill in any missing data points, and then run dPCA subsequently. Second, the number of neurons needs to be sufficiently high in order to obtain reliable estimates of the demixed components. In our datasets, we found that at least ~100 neurons were needed to achieve satisfactory demixing. The number is likely to be higher if more than three task parameters are to be demixed, as the number of interaction terms grows exponentially with the number of parameters. This trade-off between model complexity and demixing feasibility should be kept in mind when deciding how many parameters to put into the dPCA procedure. In cases when there are many task parameters of interest, dPCA is likely to be less useful than the more standard parametric single-unit approaches (such as linear regression).
As a trivial example, imagine that only N ¼ 1 neuron has been recorded; it might have strong and significant tuning to various parameters of interest, but there is no way to demix (or decode) these parameters from the recorded 'population. ' Third, even with a large number of neurons, a dataset may be non-demixable, in which case dPCA would fail. For instance, if the high-variance directions of the stimulus and the decision parts of the neural activities fully overlap, then there is no linear decoder that can demix the two parameters. Finally, dPCA components corresponding to the same parameter (e. g. successive stimulus components) are here chosen to be orthogonal, similarly to PCA. This can make successive components difficult to interpret (e. g. the second and the third stimulus components in Figure  3 ).
Materials and methods We will first explain the dPCA algorithm in the most well-behaved case of simultaneously recorded and fully balanced data. A dataset with categorical predictors is called balanced when there is the same number of data points for each combination of predictors; in our case this means that there is the same number of trials for each combination of task parameters. This scenario is unlikely in most practical applications where the experimenter often does not have full control over some of the task parameters (such as e. g. animal's decisions). Our suggestion for unbalanced datasets is to use what amounts to a 're-balancing' procedure as explained below. Finally, we will deal with the case of sequentially recorded neurons (all datasets analyzed in this manuscript fall into this category).
Mathematical notation In each of the datasets analyzed in this manuscript, trials can be labeled with two parameters: 'stimulus' and 'decision'. Note that a 'reward' label is not needed, because its value can be deduced from the other two due to the deterministic reward protocols in all tasks. In this situation, for each stimulus s (out of S) and decision d (out of Q), we have a collection of K trials with N neurons recorded in each trial. For each trial k (out of K) and neuron n (out of N) we have a recorded spike train. We denote the filtered (or binned) spike train by xðtÞ, and assume that it is sampled at T time points t. To explicitly denote all task parameters, we will write either xðt; s; d; kÞ or x tsdk for the filtered spike train of one neuron and x tsdk for the vector of filtered spike trains of all N neurons.
For the purposes of demixing neural signals in the context of our datasets, we combine some of these terms together. Indeed, demixing a time-independent pure stimulus term " x s from a stimulustime interaction term " x ts makes little sense because we expect all neural components to change with time. Hence, we group the terms as follows (without changing the notation): x tsdk ¼ " x þ " x t þ " x s þ " x ts |fflffl ffl{zfflffl ffl} " x ts þ " x d þ " x td |fflfflffl{zfflfflffl} " x td þ " x sd þ " x tsd |fflfflfflffl ffl{zfflfflfflffl ffl} " x tsd þ tsdk : Here the first term on the right-hand side is the mean firing rate, the last term is the trial-to-trial noise, and we call the other terms condition-independent term, stimulus term, decision term, and stimulus-decision interaction term. This decomposition is illustrated in Figure  8  for several exemplary neurons (we only show the decomposition of the PSTH part, leaving out the noise term).
Regularization A standard way to avoid overfitting in regression problems is to add a quadratic penalty to the cost function, which is often called ridge regression (RR). This approach can be used in reduced-rank regression as well. Specifically, we can add a ridge penalty term to the loss function L f : L f ¼ kX f À FDXk 2 þ kFDk 2 : The RR solution modifies the OLS solution from above to A RR ¼ X f X > À XX > þ I Á À1 : In turn, the reduced-rank solution can be obtained as described above: F ¼ U q and D ¼ U > q A RR where U q are the first q principal directions of A RR X. We found it convenient to define ¼ ðlkXkÞ 2 , since this makes the values of l comparable across datasets. As explained below, we used cross-validation to select the optimal value of l in each dataset.
Unbalanced data The data and variance decomposition carried out by the marginalization procedure can break down when the dataset is unbalanced, i. e. , when the number of data points (trials) differs between conditions. We illustrate this problem with a two-dimensional toy example in Figure  9. We assume two task parameters (factors), each of which can take only two possible values. The overall mean as well as the interaction term are taken to be zero, so that x ijk ¼ a i þ b j þ e ijk. Since the number of trials, K ¼ K ij , depends on the condition, the trial index runs through the values k ¼ 1. K ij. As shown in Figure  9a , all three terms on the right-hand side exhibit zero correlation between x 1 and x 2. A balanced dataset with the same number of data points in each of the four possible conditions (Figure  9b ) also has zero correlation.
However, an unbalanced dataset, as shown in Figure  9c , exhibits strong positive correlation ( ¼ 0:8). Accordingly, the covariance matrix of the full data can no longer be split into marginalized covariances. To avoid this and other related problems, we can perform a 're-balancing' procedure by reformulating dPCA in terms of PSTHs and noise covariance. In the balanced case, the dPCA loss function L f can be rewritten as the sum of two terms with one term depending on the PSTHs and another term depending on the trial-to-trial variations, L f ¼ kX f À FDXk 2 ¼ kX f À FDðX À X noise Þk 2 þ kFDX noise k 2 ; where we used the fact that X f and X À X noise are orthogonal to X noise (see Appendix A). We now define X PSTH ¼ X À X noise which is simply a matrix of the same size as X with the activity of each trial replaced by the corresponding PSTH.
Sequentially recorded data For sequentially recorded datasets, the matrix X cannot be meaningfully constructed. However, we can still work with the PSTH matrix e X that can be decomposed into marginalizations: e X ¼ P f e X f. Consequently, we can use the same formulation of the loss function as in the simultaneously recorded unbalanced case (see above). The only difference is that the noise covariance matrix is not available (noise correlations cannot be estimated when neurons are recorded in different sessions). In this manuscript we took as C noise the diagonal matrix with individual noise variances of each neuron on the diagonal. We used the re-balanced version e C noise (average noise covariance matrix across all conditions), but found that the difference between re-balanced and non-rebalanced noise covariance matrices was always minor and did not noticeably influence the dPCA solutions.
Variance calculations As all datasets analyzed in this manuscript were sequentially recorded, we always reported fractions of the PSTH variance (as opposed to the total PSTH+noise variance) explained by our components, i. e. fractions of variance explained in e X. We defined the fraction of explained variance in a standard way: R 2 ¼ k e Xk 2 À k e X À FD e Xk 2 k e Xk 2 : This formula can be used to compute the fraction of variance explained by each dPCA component (by plugging in its encoder f and decoder d); these are the numbers reported on Figures  3-6b,d  and used to order the components. The same formula can be used to compute the cumulative fraction of variance explained by the first q components (by stacking their encoders and decoders as columns and rows of F and D respectively); these are the numbers reported on Figures  3-6c.
The residual noise variance is not split uniformly across marginalizations: the fraction falling into marginalization f is proportional to the respective number of degrees of freedom, K f. This can be explicitly computed; for a centered dataset with S stimuli, Q decisions, and T time points the total number of degrees of freedom (per neuron) is SQT À 1 and is split into T À 1 for time, ST À T for stimulus, QT À T for decision, and SQT À ST À QT þ T for the stimulus-decision interaction (compare with the formulas in the Marginalization Procedure section). Accordingly, we computed the residual noise sum of squares falling into marginalization f as Q f ¼ K f SQT À 1 Q: The pie charts in Figures 3-6d show the amount of variance in each marginalization, with estimated contributions of the residual noise variance subtracted: À k e X f k 2 À Q f ÁÀ k e Xk 2 À Q Á.
Demixing indices We defined the demixing index of each component as max f fkd e X f k 2 g=kd e Xk 2. This index can range from 1/4 to 1 (since there are four marginalizations) and the closer it is to 1, the better demixed the component is. As an example, for the somatosensory working memory dataset, the average demixing index over the first 15 PCA components is 0. 76AE0. 16 (meanAESD), and over the first 15 dPCA components is 0. 98AE0. 02, which means that dPCA achieves much better demixing (p ¼ 0:0002, Mann-Whitney-Wilcoxon ranksum test). For the first 15 components of dPCA-2011  (Brendel et al. , 2011)  it was 0. 95AE0. 03, significantly less than for the current dPCA (p ¼ 0:0008). This difference may seem small, but is clearly visible in the projections by the naked eye. For comparison, the average demixing index of individual neurons in this dataset is 0.
Experimental data Brief descriptions of experimental paradigms are provided in the Results section and readers are referred to the original publications for all further details. Here we describe the selection of animals, sessions, and trials for the present manuscript. In all experiments neural recordings were obtained in multiple sessions, so most of the neurons were not recorded simultaneously. All four datasets used in this manuscript have been made available at http://crcns. org (Romo et al. , 2016; Constantinidis et al. , 2016; Feierstein et al. , 2016; Uchida et al. , 2016). 1. Somatosensory working memory task in monkeys  (Romo et al. , 1999; Brody et al. , 2003). We used data from two monkeys (code names RR14 and RR15) that were trained with the same set of vibratory frequencies, and we selected only the sessions where all six frequencies {10, 14, 18, 26, 30, 34} Hz were used for the first stimulation (other sessions were excluded).
Monkeys made few mistakes (overall error rate was 6%), and here we analyzed only correct trials. Monkey RR15 had an additional 3 s delay after the end of the second stimulation before it was cued to provide the response. Using the data from monkey RR13 (that experienced a different frequency set) led to very similar dPCA components (data not shown). 2. Visuospatial working memory task in monkeys  (Qi et al. , 2011; Meyer et al. , 2011; Qi et al. , 2012). We used the data from two monkeys (code names AD and EL) that were trained with the same spatial task. Monkeys made few mistakes (overall error rate was 8%), and here we analysed only correct trials. The first visual stimulus was presented at 9 possible spatial locations arranged in a 3Â3 grid (Figure  4a ); here we excluded all the trials where the first stimulus was presented in the centre position.
3. Olfactory discrimination task in rats  (Feierstein et al. , 2006). We used the data from all five rats (code names N1, P9, P5, T5, and W1). Some rats were trained with two distinct odors, some with four, some with six, and one rat experienced mixtures of two fixed odors in varying Figure  10. Re-stretching (time warping) procedure. We defined several alignment events (such as odour poke in, odour poke out, etc. ) and for each trial found the times t i of these events. After aligning all trials on t 1 ¼ 0 (left) we computed median times T i for all other events. Then for each trial we re-stretched the firing rate on each interval ½t i ; t iþ1 to align it with ½T i ; T iþ1 (right). After such re-stretching, all events are aligned and the trials corresponding to one condition can be averaged.
DOI: 10. 7554/eLife. 10989. 015 proportions. In all cases each odor was uniquely associated with one of the two available water ports (left/right). Following the original analysis  (Feierstein et al. , 2006) , we grouped all odors associated with the left/right reward together as a 'left/right odor'. For most rats, caproic acid and 1-hexanol (shown in Figures  5-6a ) were used as the left/right odor. We excluded from the analysis all trials that were aborted by rats before reward delivery (or before waiting 0. 2 s at the reward port for the error trials). 4. Olfactory categorization task in rats  (Kepecs et al. , 2008). We used the data from all three rats (code names N1, N48, and N49). Note that recordings from one of the rats (N1) were included in both this and previous datasets; when we excluded it from either of the datasets, the results stayed qualitatively the same (data not shown).
Selection of neurons For our analysis, we only selected neurons which had been recorded in each possible condition (combination of parameters), which avoids the missing data problems explained above. Additionally, we required that in each condition there were at least K min > 1 trials, to reduce the standard error of the mean when averaging over trials, and also for cross-validation purposes. The cutoff was set to K min ¼ 5 for both working memory datasets, and to K min ¼ 2 for both olfactory datasets (due to less neurons available). We have further excluded very few neurons with mean firing rates over 50 Hz, as such neurons can bias the variance-based analysis. Firing rates above 50 Hz were atypical in all datasets (number of excluded neurons for each dataset: 5 / 2 / 1 / 0). This exclusion had a minor positive effect on the components.
We did not apply any variance-stabilizing transformations, but if the square-root transformation was applied, the results stayed qualitatively the same (data not shown). No other pre-selection of neurons was used. This procedure left 832 neurons (230 / 602 for individual animals, order as above) in the somatosensory working memory dataset, 956 neurons (182 / 774) in the visuospatial working memory dataset, 437 neurons in the olfactory discrimination dataset (166 / 30 / 9 / 106 / 126), and 214 neurons in the olfactory categorization dataset (67 / 38 / 109).
Preprocessing of the neural data The spike trains were filtered with a Gaussian kernel (s ¼ 50 ms) and sampled at 100 Hz to produce single-trial instantaneous firing rates. In the visuospatial working memory dataset we identified the preferred location of each neuron as the location that evoked maximum mean firing rate in the 500 ms time period while the first stimulus was shown. The neural tuning was shown before to have a symmetric bell shape (Qi et al. , 2011;  Meyer et al. , 2011) , with each neuron having its own preferred location. We then re-sorted the trials (separately for each neuron) such that only five distinct stimuli were left: preferred location, 45 ˚, 90 ˚, 135 ˚, and 180 ˚away from the preferred location. In both olfactory datasets trials were self-paced. Accordingly, trials last different amounts of time, and firing rates cannot simply be averaged over trials.
We used the following time warping (restretching) procedure to equalize the length of all trials and to align several events of interest  (Figure 10)  separately in each dataset. We defined five alignment events: odor poke in, odor poke out, water poke in, reward delivery, and water poke out. First, we aligned all trials on odor poke in (T 1 ¼ 0) and computed median times of the four other events T i ; i ¼ 2. 5 (for the time of reward delivery, we took the median over all correct trials). Second, we set DT to be the minimal waiting time between water port entry and reward delivery across the whole experiment (DT ¼ 0:2 s for the olfactory discrimination task and DT ¼ 0:3 s for the olfactory categorization task). Finally, for each trial with instantaneous firing rate xðtÞ we set t i ; i ¼ 1.
Cross-validation to select regularization parameter As noted above, we renormalized the regularization parameter ¼ ðlkXkÞ 2 , and then used crossvalidation to find the optimal value of l for each dataset. To separate the data into training and testing sets, we held out one random trial for each neuron in each condition as a set of SQ test 'pseudotrials' X test (as the neurons were not recorded simultaneously, we do not have recordings of all N neurons in any actual trial). Remaining trials were averaged to form a training set of PSTHs e X train and an estimate of the noise covariance matrix e C train. Note that X test and e X train have the same dimensions. We then performed dPCA on e X train for various values of l between 10 À7 and 10 À3 (on a logarithmic grid). For each l, we selected ten components in each marginalization (i.
e. 40 components in total) to obtain F f ðlÞ and D f ðlÞ, and computed the normalized reconstruction error L CV ðlÞ on the test set (see below). We repeated this procedure ten times for different train-test splittings and averaged the resulting functions L CV ðlÞ. In all cases the average function " L CV ðlÞ had a clear minimum (Figure  11 ) that we selected as the optimal l. The values of l selected for each dataset were 2:6 Á 10 À6 / 5:8 Á 10 À6 / 5:8 Á 10 À6 / 5:8 Á 10 À6. We also performed the same procedure in each marginalization separately, but in all datasets the optimal values of l were similar across marginalizations  (Figure 11). We therefore chose to use the same value of l for all marginalizations. Interestingly, for all our datasets min l f " L CV ðlÞg was only slightly smaller than " L CV ð0Þ, so the regularization term had almost no influence.
(2012). We then used 100 shuffles to compute the distribution of classification accuracies expected by chance. On each iteration and for each neuron, we shuffled all available trials between conditions, respecting the number of trials per condition (i. e. all P sd K nsd trials were shuffled and then randomly assigned to the conditions such that all values K nsd stayed the same). Then exactly the same classification procedure as above (with 100 cross-validation iterations) was applied to the shuffled dataset to find mean classification accuracy for the first stimulus, decision, and interaction components. All 100 shuffling iterations resulted in a set of 100 time-dependent accuracies expected by chance. The time periods when actual classification accuracy exceeded all 100 shuffled decoding accuracies in at least ten consecutive time bins are marked by black lines on Figures  3 4 5 6.
Components without any periods of significant classification are not shown. See Figure  12  for classification accuracies in each dataset. The Monte Carlo computations took ~8 hr for each of the larger datasets on a 6 core 3. 2 Ghz Intel i7-3930K processor. Implementation of classical approaches (Figure  1 ) The two-way ANOVA shown in Figure  1c -e was performed as follows. The two factors were stimulus (with six levels) and decision (with two levels), the interaction term was included, and a separate ANOVA was run for the firing rate of each neuron at each time point. Significance level was set at a ¼ 0:05. Effect size was defined as partial omega squared with a sign given by the sign of the correlation coefficient between firing rate and the corresponding parameter. It can take values between À1 and 1, with 0 meaning no effect.
One particular advantage of the dPCA formulation is that it operates on the raw data, so that the decoders (and encoders) can actually be used on single trials. In turn, the visualization of the data found through dPCA also provides insights into the utility of the respective population code for the brain. In the second case we obtain U f X components from the full data, so that these components could be obtained by projecting single-trial activities. However, now there is no guarantee that these components will be demixed. For a simple counter-example, consider Figure  13 : the stimulus marginalization X s consists of three points (one for each stimulus) located roughly on a horizontal axis, and so the first principal axis of X s is roughly horizontal. It is easy to see that the projection of the full data onto this axis will be not only stimulus-, but also time-dependent.