Cross-validation in high-dimensional spaces: a lifeline for least-squares models and multi-class LDA
Matthias S. Treder











































Cardiff University Brain Research Imaging Centre (CUBRIC), Cardiff University, United Kingdom



Least-squares models such as linear regression and  Linear Discriminant Analysis (LDA) are amongst the most popular statistical learning techniques. However, since their computation time increases cubically with the number of features, they are inefficient in high-dimensional neuroimaging datasets.
Fortunately, for k-fold cross-validation, an analytical approach has been developed that yields the exact cross-validated predictions in least-squares models without explicitly training the model. Its computation time grows with the number of test samples.
Here, this approach is systematically investigated in the context of cross-validation and permutation testing. LDA is used exemplarily but results hold for all other least-squares methods. Furthermore, a non-trivial extension to multi-class LDA is formally derived. The analytical approach is evaluated using complexity calculations, simulations, and permutation testing of an EEG/MEG dataset.
Depending on the ratio between features and samples, the analytical approach is up to 10,000x faster than the standard approach (retraining the model on each training set). This allows for a fast cross-validation of least-squares models and multi-class LDA in high-dimensional data, with obvious applications in multi-dimensional datasets, Representational Similarity Analysis, and permutation testing.



MVPA classification cross-validation permutation testing LDA high-dimensional spaces RSA



















§ INTRODUCTION


Multivariate pattern analysis (MVPA) is a statistical technique in which a target variable such as a brain state or reaction time is predicted based on multivariate patterns of brain activity <cit.>. 
The spadework in MVPA is performed by regression models if the dependent variables is continuous (e.g. reaction time), or by classifiers if the dependent variable is categorical (e.g. stimulus type) <cit.>. 
Due to their simplicity, relatively low computational demands, and high interpretability, least-squares models have been popular for both regression problems (linear regression, ridge regression) and for classification problems (Linear Discriminant Analysis  <cit.>).


The increase in storage capabilities, working memory, and computational power, and the increasing availability of high-performance compute clusters paved the way for large-scale analyses of neuroimaging data. Analyses can deal with larger throughput than ever before, such as higher field strengths in fMRI and larger number of electrodes in EEG, or simply a larger amount of derived features such as time-frequency and connectivity metrics. It is worth stressing that most neuroimaging datasets have a P ≫ N shape, that is the number of features P is much larger than the number of samples N. An extreme example of this is gene expression data comprising tens of thousands of genes (features) but not more than a few hundred patients (samples) <cit.>. In cognitive neuroscience, the number of samples for an analysis is naturally capped by limits of experiment time and group size. For level 1 analyses, the number of trials is limited by the amount of time the subject can spend in the scanner. Even in fast-paced EEG/MEG experiments, it is very rare that more than 10,000 trials are collected. For level 2 analyses, sample size is equal to the number of subjects, which is typically less than a few hundred. 
Summarising, the principal challenge in large neuroimaging datasets is to efficiently cope with high-dimensional data.


Unfortunately, this is exactly the Achilles heel of least-squares methods (LSM) such as linear regression, ridge regression, and Linear Discriminant Analysis (LDA). The  computationally most expensive part in LSM is the inversion of the features × features scatter matrix. Computation time increases cubically with the number of features. It can therefore be intractable for even a few thousand features if a large number of training-testing iterations is needed, such as in permutation testing or in Representational Similarity Analysis <cit.> with many experimental conditions. This is one of the reasons that some researchers explore kernel methods such as Support Vector Machines <cit.> whose complexity grows with the number of samples rather than number of features.


Does this mean that, for all practical purposes, high-dimensional datasets are beyond reach for LSM? Fortunately, for cross-validation <cit.>, an alternative has been developed that addresses this issue. The analytical approach for LSM has the following  <cit.>useful property: instead of requiring the inversion of a feature× features scatter matrix on each training set, it instead relies on the inversion of a matrix that grows with the number of test samples. It is therefore only mildly affected by the number of features. 
For leave-one-out cross-validation, the analytical approach is well-known in the linear regression literature
<cit.>, and it has been generalised to k-fold cross-validation <cit.>. 


The aim of this study is to show that LSM can successfully meet the challenges of high-dimensional data when using the analytical approach. Because of the formal equivalence between linear regression and LDA (resp. ridge regression and regularised LDA), it suffices to focus on LDA alone. All results automatically generalise to linear regression and ridge regression. Equipped with regularisation techniques such as ridge regularisation <cit.> or shrinkage regularisation <cit.>, LDA is robust to overfitting in high-dimensional data. It often performs similarly to more sophisticated linear classifiers such as linear support vector machines (SVM) while being significantly faster to train <cit.>. 


The novel contributions in this manuscript are a detailed empirical evaluation of the analytical approach for cross-validation using simulations and complexity calculations. Furthermore, to the best of my knowledge, this is the first application in permutation testing which is a popular approach in statistical testing of classifier performance <cit.>. It is also the first time that the approach is formally extended to multi-class LDA using an optimal scoring approach <cit.>.


The manuscript is structured as follows. Firstly, cross-validation, binary LDA and a regression formulation of LDA are introduced. Then  the analytical approach is developed for cross-validation and permutation testing, and both ridge regularisation and shrinkage regularisation are considered. Finally, it is formally extended to multi-class LDA <cit.>. The computation time of the approach is then compared to the computation time of the standard approach (retraining the model for every fold) using complexity calculations, simulations, and a permutation analysis of a publicly available EEG/MEG dataset <cit.>.






§ METHOD


Matrices will denoted by bold upper case letters, for instance . Vectors are denoted as bold lower case letters, , and are assumed to be column vectors. For scalars, normal font type is used. Upper case scalars are used for specifying the dimensionality of the matrices or vectors. N is the number of samples, P the number of features or predictors, K is the number of cross-validation folds, C is the number of classes, and T is the total number of training-testing iterations in during permutation testing.







 §.§ Cross-validation




Cross-validation allows to estimate predictive performance while at the same time controlling for overfitting and making efficient use of the available samples <cit.>.

In k-fold cross-validation, the dataset is randomly partitioned into K equally sized folds. The classifier is trained on all but one of the folds, and then tested on the held out fold. This procedure is repeated until every fold served as test set once. Classification performance is then averaged across the test folds. To reduce the variance stemming from the random partitioning of data into folds, the cross-validation can be repeated several times, finally averaging across the repeats. 






 §.§ Linear Discriminant Analysis (LDA)


For two classes, LDA is equivalent to Fisher Discriminant Analysis (FDA) <cit.>. The multi-class case is considered further below. Geometrically speaking, LDA seeks a projection  from feature space to a 1-dimensional subspace such that the projected class means are maximally separated while at the same time the projected variance within classes is minimised <cit.>. Using the LDA derivation in Duda & Hart <cit.> this can be formalised as:


    J() = ^⊤_b /^⊤_w


where _b∈^P× P is the between-classes scatter matrix measuring the distance between the classes and _w∈^P× P is the within-class scatter matrix measuring the spread within each class. These quantities are defined as


    2_b =    ∑_l ∈{1,2}N_l (l -) (l - )^⊤         (between-classes scatter)
    _w =    ∑_l ∈{1,2}∑_i∈𝒞_l (_i - l)(_i - l)^⊤         (within-class scatter)


where 𝒞_l is an index set representing the samples in class l, N_l is the number of samples in class l, and the means are given by


    2l =    1/N_l∑_i∈𝒞_l_i         (class mean)
     =    1/N∑_i∈{1,2,...,N}_i         (sample mean)


Obviously, the sample mean and the two class means are related by N  = N_11 + N_22.
In the two-classes case, _b further simplifies to


    _b = N_1 N_2/N (1-2) (1-2)^⊤  (between-classes scatter)


Setting λ = J() one arrives at the generalised eigenvalue problem _b  = λ _w. For a binary classification problem and a positive definite within-class scatter matrix, the eigenvector corresponding to the largest eigenvalue is proportional to


    = _w^-1 (1 - 2)  (weight vector)


This is proved in lem:evproblem in app:proofs. It is expedient to define the bias as the center between the projected class means because this prevents the classifier from being biased towards one of the classes if the number of training samples per class is not equal:


    b_LDA = -^⊤ (1 - 2)/2  (bias term)


The classifier output ŷ for a new sample  is calculated as ŷ := ^⊤ + b. This quantity is the signed distance to the hyperplane, more generally known as decision value. It is these decision values that are subject to the analytical approach developed below. The class labels can be derived from the sign of the decision value, with class "+1" for ŷ≥ 0 and class "-1" for ŷ<0 <cit.>. 






 §.§ Binary LDA as a least-squares problem






There are several equivalent formulations of LDA.
For two classes, LDA is formally equivalent to LCMV beamforming <cit.>. Furthermore, as illustrated in fig:LDA_as_regression, binary LDA can be cast as a least-squares regression problem <cit.>. Let ∈^N× P be the data matrix containing samples as rows and features as columns. Then  ∈^N× (P+1) is the augmented data matrix  obtained by adding a column of 1's (for the bias term) and β is the regression weights vector absorbing both  and b 


    = [, 1_N]∈^N× (P+1),   β = [     ; b_LR ]∈^P+1


where _N is a vector of N ones. The bias term is denoted as b_LR since it is generally different from the LDA bias term b_LDA. The class labels are collected in a response vector ∈^N that uses numerical codes (e.g. +1 and -1) for the class labels. The standard regression problem using the full dataset can then be formulated as


    Β̂ = Βarg min || Β - ||_2^2


with the solution given by Β̂ = (^⊤)^-1 ^⊤. 




Following the derivation in app:app-regression-binary-fda, and assuming that the class labels are coded as +1 and -1, it can be shown that Β̂ consists of the two components


    ∝    _w^-1 (_1 - _2)
    
    b_LR =    N_1 - N_2/N - ^⊤


In other words, the solution for  using linear regression is proportional to the LDA solution given in eq:w. Since scaling does not affect classification performance, one can say that the solutions are identical. Unless the classes have equal proportions of samples (N_1=N_2), the bias term b_LR differs from the common choice presented in eq:b. However, class proportions and the exact numerical coding of the classes in  do not affect the direction of . This is shown in app:app-regression-binary-fda.























 §.§ An analytical approach to cross-validation for least-squares methods




In this section, the analytical approach to cross-validation is introduced for binary LDA. For linear regression and ridge regression, the approach is identical, with  simply being continuous a response variable instead of a vector of class labels. The approach has been introduced before for leave-one-out cross-validation <cit.> and k-fold cross-validation <cit.> but without the detailed derivation provided here. 

Let Tr ⊂{1,2,…,N} be the indices of the training samples, and Te ⊂{1,2,…,N} be the indices of the test samples.

Let ∈{-1,+1}^N be the vector of class labels. In the regression framework, these class labels serve as target decision values. The decision values obtained from the classifier trained on the full dataset are denoted as . The cross-validated decision values obtained from a classifier trained on only the training set and then tested on the independent test set are denoted as . _ resp. _, and _ resp. _ refer to the submatrix or subvector corresponding the training resp. test samples. To ease reading of the formulas, some important quantities used in the derivations are depicted in fig:symbol_explanation.



  §.§.§ Basic idea


In the regression framework, training the classifier is equivalent to calculating the vector of regression weights  from the training data


    = (_^⊤ _)^-1 _^⊤ _  (model based on training data)


In k-fold cross-validation, this process is repeated for each of the K training folds. However, as will be shown next, it suffices to train only one model using all available data


    = (^⊤ )^-1 ^⊤   (model based on all data)


and then obtain the cross-validated decision values directly via an analytical approach.



  §.§.§ Hat matrix


The hat matrix ∈̋^N× N is defined as


    =̋ (^⊤)^-1 ^⊤  (hat matrix)


It is a quantity well-known in linear regression <cit.>. Its name stems from the fact that it "puts the hat" onto the response vector  by mapping the true responses onto the predicted responses =. The submatrix  is obtained from $̋ by selecting only the rows and columns that correspond to test samples.is obtained from$̋ by selecting the rows corresponding to training samples and the columns corresponding to test samples.
As will be seen below, the hat matrix arises naturally during updating.



  §.§.§ Updating _^⊤ _


A formula that will prove useful later on is the product _^⊤ _ which can be obtained as follows


    _^⊤ _ = ^⊤ - ^⊤ _.





  §.§.§ Updating the inverse scatter matrix


Similarly, the scatter matrix on the training data can be obtained from the full scatter matrix by removing the scatter corresponding to the test samples, 


    _^⊤ _ = ^⊤  - ^⊤ .


Suppose that ^⊤ and its inverse, denoted as := (^⊤ )^-1 have already been calculated and let  be the identity matrix. The matrix inversion lemma (a.k.a. Sherman-Morrison-Woodbury formula) can be used to update the inverse scatter matrix on the training data as


    (_^⊤ _)^-1 =    (^⊤  - ^⊤ )^-1
    
    =     +  ^⊤ ( -   ^⊤)^-1 
    
    =     +  ^⊤ ( - )^-1


This solution circumvents the explicit inversion of the scatter matrix, but it still involves a number of matrix multiplications. It therefore only serves as an intermediate result.



  §.§.§ Updating the weights


To calculate the weights on the training data, eq:update-Xy and eq:update-inverse  can be plugged into eq:Beta-train. This yields


    =    ( +  ^⊤ ( - )^-1 ) (^⊤ - _^⊤ _)
    
    =     - ( ^⊤ ( - )^-1) ( [-] _ -   +  _)
    
    =     - ( ^⊤ ( - )^-1) ( _ -  )
    
    =     - ( ^⊤ ( - )^-1) ( _ - _)


where _  := _ - _ is the estimation error on the test samples between the correct and the predicted decision values using a model trained on the full dataset.







As will be seen next,  does not need to be calculated explicitly. 



  §.§.§ Updating the decision values


The goal is to derive the cross-validated decision values for the test samples denoted as _. As an intermediate step, the corresponding estimation error is calculated first.


    _= _ - _  (cross-validated estimation error)


Inserting eq:update-beta then leads to the desired analytical approach


    _ =    _ - _ 
    
    =    _ - __=_ +  ( - )^-1 _
    
    =    ( -  + ) ( - )^-1 _
    
    =    ( - )^-1 _    (analytical approach)


It is now easy to obtain the cross-validated decision values on the test set by simply solving eq:ecvte for _. Finally, these decision values can be used to calculate classification accuracy, AUC, or any other desired metric of classification performance.



 §.§ Adjusting the bias term


The bias term resulting from the regression approach does not generally coincide with the bias term used in LDA. However, the bias terms to coincide if N_1 = N_2. Hence, for unbalanced data, undersampling of the majority class or oversampling of the minority class is a remedy. Alternatively, if area under the ROC curve (AUC) is used as classifier performance metric, the bias term is irrelevant. 

If it is not possible to use one of these approaches, the bias term needs to be adjusted. To this end, the class means and the sample mean on the training need to be calculated and projected onto  (see definition of b_LDA and b_LR). Fortunately, it is not required to explicitly calculate . Instead, one can determine the decision values of the cross-validated model on the training set, _, and then calculate b_LR and b_LDA directly. This is achieved by applying eq:update-rule to the training data:


    _ =    _ + _̋,  ( - )^-1 _
    _ =    _ - _


Finally, the operation _←_ - b_LR + b_LDA adjusts the bias.



  §.§.§ Summary: analytical approach


It has been shown that in k-fold cross-validation, it suffices to train just one regression model on the whole dataset. By evaluating eq:update-rule the decision values for each of the folds are obtained directly without explicitly calculating any of the K models.





 §.§ Regularisation


Since neuroimaging data is often low-dimensional, or the number of samples is smaller than the number of features, the within-class scatter matrix tends to be ill-conditioned. Two similar regularisation approaches have been explored in the literature. They are equivalent in that they define the same family of classifiers (up to scaling of ).



  §.§.§ Ridge regularisation


A multiple of the identity matrix is added to the within-class scatter matrix <cit.>. The regularised within-class scatter matrix is  _w + λ, 






where λ∈[0,∞] is the regularisation term and  is the identity matrix. λ = 0 yields the ordinary, unregularised solution. The ridge solution for  is then given by


    = (_w + λ)^-1 (1 - 2)


In app:ridge it is proven that, in the regression framework, the corresponding solution is given by


    Β̂ = (^⊤ + λ_0)^-1 ^⊤


where _0∈^(P+1)×(P+1) is a diagonal matrix that is identical to the identity matrix except for the last element which is 0 instead of 1. This construction assures that the bias term corresponding to the last entry is not subjected to regularisation.

Analogous to eq:update-scatter, the regularised scatter matrix for the training data can be obtained as an update on the full model:


    _^⊤ _ + λ_0 = ^⊤  + λ_0_full scatter - ^⊤ _update


After redefining := (^⊤  + λ_0)^-1 and correspondingly including the regularisation term in the hat matrix =̋ (^⊤ + λ_0)^-1 ^⊤, the analytical approach is identical to eq:update-rule.



  §.§.§ Shrinkage regularisation


The within-class scatter matrix is replaced by a convex combination of the empirical covariance matrix and a scaled identity matrix, (1-λ) _w + λν, where ν = trace(_w)/P is a scaling parameter that equalises the traces of _w and ν, and λ∈[0,1] <cit.>. Unfortunately, shrinkage regularisation does not allow for simple low-rank updates as before. This can be seen when one inspects an update of the regularised scatter matrix


    (1-λ)_^⊤_ + λν__0 = (1-λ)^⊤ + λν_0_full scatter - ((1-λ)^⊤ + λ(ν-ν_)_0)_update


where ν is the scaling calculated on the full dataset and ν_ is the scaling calculated on the training data. The problem is that it is necessary to update the regularisation term as well. This is caused by the scaling factor ν_, which changes for each training set, thereby changing the amount of regularisation. This turns a low-rank update into a full rank update, precluding significant performance gains by updating.

For this reason, it is recommended to resort to ridge regularisation. If a researcher is used to work with shrinkage, the following simple relation can be used to transform a given shrinkage parameter λ_shrink into a corresponding ridge parameter λ_ridge. Given a fixed value for λ_shrink, the goal is to find a λ_ridge such that the regularised scatter matrices are proportional:


    (1-λ_shrink) ^⊤  + λ_shrink ν _0  !∝  ^⊤  + λ_ridge _0


Obviously, this relation holds when the ridge parameter is defined as


    λ_ridge = λ_shrink/1-λ_shrink ν







 §.§ Using the analytical approach for permutation testing


The hat matrix $̋ is invariant under class label permutations because it depends on the features alone. Consequently, it does not need to be recalculated when the class labels are permuted. Let the permuted class labels be denoted as. Ifandare adjusted accordingly


    ←    
    ←


the formulas in the previous section directly apply. They are compiled in alg:alg1.






 §.§ Multi-class LDA


Multi-class LDA is the generalisation of binary LDA to more than two classes. Like binary LDA, it involves a projection step and a thresholding step. In the projection step, the data is mapped onto a(C-1)-dimensional subspace, whereCis the number of classes. In the second step, a new sample is assigned to the class with the closest class centroid. LDA thus acts as a prototype classifier.
The scatter matrices are calculated as before, but now information is pooled across all classes.


    2_b =    ∑_j ∈{1,2,...,C}n_j (j -) (j - )^⊤         (between-classes scatter)
    _w =    ∑_j ∈{1,2,...,C}∑_i∈𝒞_j (_i - j)(_i - j)^⊤         (within-class scatter)


Assuming that there are more features than classes, the between-classes scatter matrix_bhas rankC-1. Consequently, there are multiple non-trivial solutions that again can be obtained via the generalised eigenvalue problem


    _b  = _w Λ


whereΛis a diagonal matrix of eigenvalues. A set of discriminant coordinates is obtained corresponding to non-zero eigenvalues of the eigenvalue problem. They are collected in a matrix∈^P×(C-1)and scaled such that^⊤_w= <cit.>. As for binary LDA, ridge regularisation can be applied to the within-class scatter matrix by replacing_wby_w+λ<cit.>.






 §.§ Multi-class LDA in a regression framework


Unfortunately, multi-class LDA is not equivalent to multivariate linear regression using a class indicator matrix as response matrix <cit.>. Nevertheless, there is a close relationship between both approaches  <cit.>.
A useful characterisation is given in <cit.>. They show that multi-class LDA is equivalent to optimal scoring (OS) wherein regression is performed using a response vector with optimal numerical scores for each class. Finding the optimal scores is an optimisation problem that is solved jointly with the regression problem.
Let∈^N×Cbe the class indicator matrix whose (i,j)-th element is defined as


    _ij =
       
         1    if sample i belongs to class j
    
         0    otherwise


Letθ∈^Cbe the vector containing the optimal scores. Then the response vector of optimal scores can be written asθ, and the optimal scoring problem is given by


    Β,θarg min || Β - θ||_2^2   (optimal scoring)


whereΒandθare jointly optimised. The additional constraintN^-1||θ||^2 = 1avoids trivial solutions.
<cit.> show that this optimisation problem can be broken up into two successive steps.
Step 1: A multivariate regression is performed on the class indicator matrix = arg min ||  - ||_F^2, 
where||·||_Fis the Frobenius norm. The result is a matrix of regression weights∈^(P+1)×C, where each column of regression weights corresponds to the respective column of. This yields the matrix of regression fits= . 
Step 2: The optimal score vector is found via an eigendecomposition of^⊤. Let∈^C×(C-1)be the eigenvectors of this decomposition, also called optimal scores, where the column corresponding to the trivial eigenvalue 0 (ifis centered) or 1 (ifis not centered) has been removed. Letα_1^2,α_2^2,...,α_C-1^2be the corresponding eigenvalues. Letbe the submatrix ofwith the last row (bias term) omitted. Then the columns ofpoint in the same directions as the discriminant coordinates obtained in multi-class LDA but their scaling differs. To scale the discriminant coordinates, they are right-multiplied with the diagonal matrix


    𝐃 = √(N)^-1[     √(α_1^2(1-α_1^2))                     0                     …                     0;                     0     √(α_2^2(1-α_2^2))                     …                     0;                     ⋮                     ⋮                     ⋱                     0;                     0                     0                     … √(α_C-1^2(1-α_C-1^2)) ]^-1


Note that the normalisation√(N)^-1does not appear in the original definition of the scaling matrix (<cit.>, p. 83) but is necessary here because the multi-class LDA has been calculated using the within-class scatter matrix. In contrast, Hastie et al. use the covariance matrix which differs by a scaling factor ofN. As main result of their derivation, the relationship between the  discriminant coordinatesin eq:LDA-eigenvalue-multiclass and the optimal scoring results is given by


    =





 §.§ The analytical approach for multi-class LDA


How can these findings be used to develop an analytical approach for multi-class LDA?


Starting from eq:LDA-as-OS, a dot is used to indicate that the matrices have been estimated using the training data. Left-multiplication with the test data then yields


    _ Ẇ =    _ 
    ⇔ _ =    _


where_is used to denote the desired discriminant scores for the test data (obtained in step 2). This notation is necessary to differentiate them from the cross-validated regression fits(obtained in step 1).  After calculating= ,_and_can be obtained by applying eq:update-rule and eq:fixb using the matrix of estimation errorsÊ = - .andare  then obtained via the eigenanalysis(_^⊤  _/ N_)on the training data.
Note that in practice, the augmented data matrixand the regression weightscan be used. The classification results are equivalent since the distance of a sample to the class centroids is unaffected by the constant shift incurred by the bias term.


Concluding, an analytical approach for step 1 of OS has been developed. There is no straightforward way to update the eigenvalue decomposition in step 2. However, the eigenanalysed matrix is of dimensionsC×C, so for most practical applications the computational costs are negligible. alg:alg-multiclass compiles these results.






 §.§ Computational complexity of the analytical approach


The asymptotic computational complexity for classifier validation using cross-validation is quantified in terms of floating point operations. The analytical approach developed in this paper is compared to the standard approach wherein a classifier is trained from scratch on every training set. Regularisation is not considered separately since the addition of the regularisation term inflicts negligible costs and is the same in both algorithms. For simplicity, the complexity is given in terms of the standard textbook algorithms. Speed-ups can of course be achieved using more sophisticated algorithms. The complexity calculations are summarised in tab:complexity.



  §.§.§ Binary LDA


For training a single binary LDA classifier, two class means need to be calculated which involves adding up training samples and features and dividing two times, leading to a complexity of NP, wherePis the number of features. Calculating the within-class scatter matrix requiresN(P+P^2)steps, where thePis for subtracting the class mean and theP^2is for calculating the outer vector product (NP^2). Instead of then calculating the inverse of the within-class scatter matrix, one can solve the system of linear equations_w = _1 - _2(P^3). Calculation of the biasb_LDArequires P. Taken together, the complexity for training a single classifier amounts to NP^2 + P^3. This process is repeatedKtimes,Kbeing the number of folds.

This yields an overall complexity of KNP^2 + KP^3 since the lower-order terms can be ignored for asymptotic complexity.

For the analytical approach based on the regression approach, the hat matrix needs to be calculated initially at a complexity of N^2P + NP^2 + P^3. Then eq:update-rule needs to be evaluated for each training iteration  (KN^3). If the bias term needs to be corrected, operations at KN^2 are required. This yields an asymptotic complexity of KN^3 for the analytical approach.



  §.§.§ Multi-class LDA


The asymptotic complexity for training a single multi-class LDA classifier is provided first. Calculating means for each of theCclasses, involves adding up training samples and features and dividingCtimes (NP + CP).
Calculating the within-class scatter matrix is equal to the binary case (NP^2). Calculating the between-classes scatter matrix involves calculatingCouter vector products CP^2. The generalised eigenvalue decomposition has an overall complexity of P^3.  Repeating this processKtimes yields an overall complexity of KNP^2 + KCP^2 + KP^3.

For the analytical approach based on the optimal scoring approach, the hat matrix needs to be calculated initially at a complexity of N^2P + NP^2 + P^3. Obtaining the cross-validated regression fits_and_involves a complexity of KN^3 C each. This is followed by the calculation and eigendecomposition of_^⊤  _(KC^2 N + KN^3 C). Finally the discriminant scores are calculated (KC^2 N). This yields an asymptotic complexity of KN^3 C.






 §.§ Simulations


To vet the analytical approach its efficacy is compared to the standard approach using simulated data. The data is created as follows: Each class centroid is randomly placed on the surface of a unit hypersphere in feature space. A common covariance matrix is randomly sampled from a Wishart distribution. Samples are then created by randomly sampling from a multivariate normal distribution parameterised by the corresponding class centroid and the common covariance matrix.


The number of features was varied from 10 to 1000 in 40 logarithmic steps. For binary LDA, cross-validation was performed using 5 folds, 10 folds, 20 folds, and leave-one-out. Simulations were run separately for 100 and 1000 samples. Simulations were run with 10-fold cross-validation and 100, 1000, or 10,000 permutations. The number of samples and the number of features were set to either 100 or 1000. For every combination of parameters,  the simulation was repeated 20 times.


For multi-class LDA, 10-fold cross-validation was used with data being split into 5 classes or 10 classes with equal class proportions. For cross-validation, the number of samples was either 100 or 1000. For permutations, the number of features was fixed to 100 or 1000. The number of permutations was limited to 10 or 100 to keep overall  computation time tractable. For every combination of parameters, the simulation was repeated 20 times for cross-validation. For permutations, it was repeated 10 times.

Both binary and multi-class LDA update rules were compared to the vanilla approach wherein the classifier is trained on each training set and then applied to the test set. Different datasets and different folds were randomly created for each choice of parameters. However, for each of the two methods (analytical approach vs classical approach) the random seed was reset to assure equal data and equal folds.
All analyses were performed in MATLAB (Natick, USA). The  and  functions were used to measure the total computation time for cross-validation and permutation testing iterations. As target measure, relative effiency was computed, defined as


    Relative efficiency = log_10 time(standard approach)/time(analytical approach)


This quantity has a simple interpretation in terms of orders of magnitude of improvement in computation time of the analytical approach over the standard method. For instance, a value of 0 means that both methods are at parity. A value of 1 means that the analytical approach is 10 times faster than the standard method, a value of 2 means that the analytical approach is10^2 = 100times faster, and so on. Simulations were run on a Thinkpad X1 Carbon with 16 GB of RAM and an Intel Core i7-6600U CPU @ 2.60GHz×4 processor.








 §.§ EEG data


The analytical approach was applied to a publicly available multi-modal dataset of participants watching greyscale images of faces and scrambled faces <cit.>. The 16 EEG/MEG datasets with a total of 380 EEG/MEG channels were read into MATLAB using FieldTrip <cit.>. Epochs were created from -0.5 s to 1 s relative to image onset, and the pre-stimulus interval was used for baseline correction. Finally, data were downsampled to 200 Hz. The type of stimulus (face vs scrambled) was used as class label for binary LDA. For multi-class LDA, the face stimuli were further split in order to create 3 classes.
The total number of  trials varied across subjects, with 787 trials on average.

For both binary and multi-class LDA, two different analyses were conducted. In the first case, classification was performed separately for every time point across time interval  ranging from -0.5 s to 1 s. At each time point, 100 permutations were conducted with shuffled class labels and a 10-fold cross-validation in each permutation. The amplitudes in each channel and were used as features (380 features). In the second case, the post-stimulus interval was divided into successive, non-overlapping windows. The amplitudes in each channel were averaged within these windows and then all averaged amplitudes were concatenated to a single feature vector. For binary LDA, 100 ms windows were used (10 * 380 = 3800features). For multi-class LDA, 200 ms windows were used (5 * 380 = 1900features).

Analyses were run on the high-performance cluster at the Cardiff University Brain Research Imaging Centre (CUBRIC). Each compute node consists of 12 cores with 192 GB RAM and an Intel(R) Xeon(R) X5660 CPU running at 2.80GHz.




 §.§ Software


MATLAB implementations of the analytical approach and the scripts reproducing the results are publicly available on GitHub (). Note that parts of the code require the MVPA-Light toolbox () to run.






§ RESULTS


 §.§ Simulations
Binary LDA. A three-way analysis of variance (ANOVA) was run on the cross-validation analysis (fig:simulation_results, top left). A continuous variable (features) and two categorical variables, samples N (100 or 1000) and folds (5, 10, 20, leave-one-out), were used as predictors, and relative efficiency was used as dependent variable.  There were significant main effects of features (F = 32051.69; p < .001), N (F = 1316.19; p < .001), and folds (F = 3119.03; p < .001). Furthermore, the effect of features increased with folds (features×N,F = 806.49; p < .001), and there was an N×folds interaction (F = 812.7; p < .001). Furthermore, there was a three-way interaction N×folds×features (F = 37.16; p < .001).

A separate three-way ANOVA was performed on the permutations data (fig:simulation_results, top right) using N, permutations, and features as predictors. There were significant main effects for N (F = 6899.92; p < .001), permutations (F = 4.35; p = .014), and features (F = 111506.59; p < .001). Significant interactions were   N×permutations  (F = 26.52; p < .001) and N×features (F = 273.66; p < .001), illustrating that the effects of permutations and features were larger for N=1000 than for N=100. Other interactions were not significant (permutations×features,p = .58; N×permutations×features,p  = .08).
Multi-class LDA. A three-way analysis of variance (ANOVA) was run on the cross-validation analysis (fig:simulation_results, bottom left). A continuous variable (features) and two categorical variables, samples N (100 or 1000) and classes (5, 10), were used as predictors, and relative efficiency was used as dependent variable.
There were significant main effects for N (F = 1023.97; p < .001), features (F = 38270.22; p < .001) but not for classes (p = .15). There was a significant 
 features×N interaction (F = 125.74; p < .001) signifying a smaller effect of features for larger N. The other interactions were not significant (features×classes,p = .1; N×classes,p = .86, features×N×classes,p =.462). 

A separate three-way ANOVA was performed on the permutations data (fig:simulation_results, bottom right) using N, permutations, and features as predictors. There were significant main effects of N (F = 366.2; p < .001), permutations (F = 27.4; p < .001), and features (F = 16970.31; p < .001).
Again, there was a significant N×features interaction (F =  24.93; p < .001) signifying a smaller effect of features for larger N. The other interactions were not significant (N×permutations,p = .13; permutations×features,p = .35, N×permutations×features,p =.06).



 §.§ EEG/MEG data


Results are depicted in fig:WakemanHenson. The data for binary LDA and multi-class LDA were combined into a single two-way ANOVA model. Features (small = 380, large = 3800/1900) and type of classifier (binary LDA, multi-class LDA) were used as predictors. There were significant main effects of features (F = 826.04; p < .001) and classifier (F = 388.2; p < .001). Moreover, there was a significant features×classifier interaction (F = 6.01; p = .017) signifying a larger effect of features for multi-class LDA than for binary LDA.







§ DISCUSSION


Due to its robustness and competitive speed, regularised LDA is an excellent candidate for classification problems involving many training-testing iterations. The analytical approach explored in this paper boosts the performance of least-squares models and multi-class LDA, particularly for high-dimensional data.


The simulations revealed a persistent speed gain using the analytical approach to cross-validation as opposed to the standard approach wherein the classifier is retrained on every training fold. Relative efficiency increases notably with the number of features. Furthermore, it increases with the number of cross-validation folds, and it decreases when the number of samples increases. Cross-validation was significantly faster, by up to 3 orders of magnitude (1000x faster) for binary LDA, and close to 4 orders of magnitude (10,000x faster) for multi-class LDA.


In line with this, the analytical approach consistently outperformed the standard approach in the EEG/MEG analysis. This was particularly prevalent in multi-class LDA, where for 1900 features, the analytical approach was between 1000x and 10,000x faster than the standard approach.


Why is relative efficiency higher for multi-class LDA than for binary LDA? A possible explanation is that multi-class LDA is more involved computationally since a generalised eigenvalue problem needs to be solved, whereas binary LDA requires only matrix inversion. Crucially, the analytical approach requires only a single matrix inversion in both  the binary and the multi-class case. This warrants a larger computational benefit for multi-class LDA.





 §.§ Is it just a trade-off between samples and features?

It is worth noting that the analytical approach does not simply trade in the number of samples (N) for the number of features (P). If both quantities are equal, e.g. N = P = 1000, binary LDA is about 10x faster than the standard approach for 10-fold cross-validation and about 100x faster for leave-one-out. For multi-class LDA and 10-fold cross-validation, relative efficiency is close to 3 (almost 1000x faster).


The complexity calculations of the analytical approach grows cubically with the number of test samples whereas the standard method grows cubically with the number of features.
Consequently, the standard method and the analytical approach are at parity when
the number of test samples is roughly equal to the number of features, i.e. N / K ≈P. One can deduce the rule of thumb that it is beneficial to use the analytical approach in cross-validation as soon asP > N / K. The approach becomes more efficient as K increases with the upper limit being leave-one-out (K = N).




 §.§ What is the practical use of the analytical approach?


The analytical approach yields a significant increase in speed, but is this practically relevant for typical neuroimaging analyses? There are a number of scenarios in modern neuroimaging analyses wherein a large number of training-testing iterations is performed and hence the approach developed here can be useful:



  * Multi-dimensional data. Sometimes, statistical analysis is repeatedly performed along multiple data dimensions. For instance, in time-frequency data, a classifier may be validated for every combination of time point and frequency. In time generalisation, a classifier is trained and tested at every combination of time points in a trial. In searchlight analysis <cit.>, a classifier is validated on a local neighbourhood centered on a voxel, and this operation is repeated for all voxels.

  * Condition-rich designs. Some experimental designs, often used in the context of Representational Similarity Analysis (RSA) <cit.>, feature a large number of stimulus conditions. To build the Representational Dissimilarity Matrix, distances between each pair of conditions are required. Hence with C conditions, C(C-1)/2 cross-validations are required for every subject. Initially RSA was based on simple Pearson correlation between samples, but more recent work has increasingly focused on classifier-based approaches, including LDA classification accuracy and LDA-related measures such as Linear Discriminant Contrast (LDC) <cit.>. 

  * Permutation testing. For permutation testing, a classification regime needs to be repeated thousands of times. For instance, <cit.> developed a cluster test for fMRI data  that involved repeated 100 classification analyses with permuted class labels for each searchlight position and every subject. The results were passed on to the second level to perform group inference. Similarly, in <cit.>, permutations at the subject level are computed for deriving a minimum-statistic used in group inference.



 §.§ LDA vs. other least-squares approaches


Although the analyses presented in this paper focus on LDA, all results readily extend to other least-squares methods such as linear regression and ridge regression. If the vector of class labels is replaced by a vector of continuous responses, then all equations and results apply equally. Furthermore, since multi-class LDA is closely related to Canonical Correlation Analysis (CCA) <cit.>, speed-ups for CCA might be possible using a similar approach. 




 §.§ Analytical approach vs. kernel methods


Kernel-based methods such as Support Vector Machines <cit.> are based on a samples×samples kernel matrix𝐊that can be thought of as representing pair-wise similarities between samples according to some non-linear similarity measure. The calculation of the kernel matrix is affected by the number of features, but once the matrix it is available, optimisation algorithms such as Dual Coordinate Descent for SVM <cit.> directly operate on the samples dimension rather than the feature dimension. For linear SVM, the kernel matrix comprises standard dot products,𝐊_ij = _i^⊤_j. There is a close relationship between the linear kernel and the hat matrix which consists of the entries_̋ij = _i^⊤(^⊤+ λ)^-1 _j. Forλ>0, this quantity is positive-definite and hence a valid dot product. In other words, the hat matrix is simply a linear kernel whereby the samples have been pre-whitened with respect to the regularised scatter matrix. If the covariance of the samples is normalised and spherical, we have=̋ 𝐊.

However, kernel methods such as SVM require iterative optimisation algorithms as well as optimisation of hyperparameters. The analytical formula for least-squares methods and multi-class LDA hence yields a computational advantage in many cases.




 §.§ What about big data?


Due to increasing levels of data sharing and large-scale studies, cognitive neuroscience is on the verge of becoming a big data science <cit.>. In a big data setting, both the number of samples and the number of features is extremely large. This poses a challenge to all statistical learning approaches. For least-squares models, this challenge is admittedly not resolved with the present contribution. However, the following measures can be used to cope with either too many samples or too many features.


  * Too many samples. The principal problem is that for a very large number of samples (e.g. >100,000) it might be impossible to store the hat matrix in memory. Since a kernel matrix has the same size as the hat matrix, a similar problem occurs in the optimisation of SVMs. In SVM, on-the-fly calculation of the required kernel matrix entries has been proposed as a solution <cit.>. If the number of features is small enough, the  matrix ^⊤ can be stored in memory and submatrices of the hat matrix can be calculated on the fly. Furthermore, the submatrices - _̋ that need to be inverted   are roughly of size N/K. Hence, one can always find a K large enough such that these matrices are small enough to be invertible efficiently.

  * Too many features. If the number of features is too large, it is impossible to store the scatter matrix ^⊤ in memory. Random projections can offer a solution to this problem. There is evidence that if ∈^N× P is multiplied by a sparse matrix 𝐀∈^P× Q with Q ≪ P, the covariance structure of the original data is approximately preserved in the smaller, sparsified matrix 𝐀∈^N× Q<cit.>. This matrix can then be used instead of the scatter matrix.


An alternative approach that deals with both issues simultaneously is ensemble learning <cit.>, wherein a large number of statistical models called weak learners is trained in parallel. Each model uses a subset of features and a subset of samples. If these subsets are small enough, even large datasets can be digested by the ensemble. Furthermore, since each weak learner is trained independently of the others, ensemble learning can be efficiencly parallelised on compute clusters.



 §.§ Conclusion

For least-squares methods and multi-class LDA, an analytical approach to cross-validation allows for an increase of computation speed up to several orders of magnitude. The analytical approach  enables least-squares methods and multi-class LDA to be used in high-dimensional feature spaces, particularly in theP ≫Nsetting (many features, few samples) often encountered in neuroimaging data. Target applications in modern neuroimaging studies include multi-dimensional datasets, Representational Similarity Analysis, and permutation testing.














§ RELATIONSHIP BETWEEN LINEAR REGRESSION AND BINARY LDA


The regression problem eq:LDA-regression-problem leads to the normal equations


    ^⊤ Β =


Recall thatis the augmented data matrix consisting of the original dataand a column of 1's. Without loss of generality, one can assume that the samples in the data matrixhave been arranged as= [_1; _2]such that samples corresponding to class 1 come first and samples corresponding to class 2 come last. The response vector∈^Ncontains the numerical codes for the class labels. Class 1 is represented byz_1∈, class 1 is represented byz_2∈, class 2 is represented byz_1≠z_2. Accordingly,consists ofN_1timesz_1followed byN_2timesz_2. Plugging this into eq:normal-equations-Xa yields


    [   _1^⊤   _2^⊤; _N_1^⊤ _N_2^⊤ ][   _1 _N_1;   _2 _N_2 ][  ; b ] =
    [   _1^⊤   _2^⊤; _N_1^⊤ _N_2^⊤ ][ z_1_N_1; z_2_N_2 ]


Multiplying the matrices and using^⊤= _w + N_1_1_1^⊤+ N_2_2_2^⊤one obtains


    [ _w + N_1_1_1^⊤ + N_2_2_2^⊤                          N;                       N ^⊤                          N ][  ; b ] =
    [ N_1 z_1_1 + N_2 z_2_2;     N_1 z_1 + N_2 z_2 ]


with_1, _2, andas defined in eq:means. Solving the last row of the equation for b yieldsb = N_1 z_1/ N + N_2 z_2/ N - ^⊤. Plugging this into the first equation in eq:normal-equations2 yields


    (_w + N_1_1^2 + N_2_2^2 - N ^2)  = N_1_1 -N_2_2 - (N_1 z_1 + N_2 z_2)


where^2is short for^⊤. UsingN= N_1_1 + N_2_2and the relationN_1 - N_1^2/N = (N_1 N_2)/None obtains


    (_w + N_1 N_2/N (_1-_2)(_1-_2)^⊤__b)  = 2N_1 N_2  (z_1-z_2)/N(_1 - _2)


The vector_b is a multiple of(_1-_2), hence there existsα∈such that


    _b  = (2N_1 N_2 (z_1-z_2)/N - α) (_1 - _2)


Inserting eq:normal-equations-alpha
 in eq:normal-equations4 and assuming that_wis regular yields


    _w = α (_1 - _2) ⇔ = α _w^-1(_1 - _2)


This proves thatin the linear regression approach is (up to scaling) identical to the LDA solution. Furthermore, the exact numerical codesz_1andz_2for the classes determineband the scaling of, but they do not affect the direction of.






§ RIDGE REGULARISATION FOR BINARY LDA


In this section, the correspondence between the regularised LDA in eq:w-ridge and ridge regression solution in eq:LDA-ridge-regression-solution is established. To simplify the math, it is assumed that in, class 1 is coded as+1and class 2 is coded as-1. The assertion is that regularised LDA can be cast in a least-squares framework using the normal equations


    (^⊤ + λ_0) Β =


where_0is defined like in eq:LDA-ridge-regression-solution. Following the derivation in the previous section, one arrives at


    (
    [   _1^⊤   _2^⊤; _N_1^⊤ _N_2^⊤ ][   _1 _N_1;   _2 _N_2 ] + λ_0
    )
    [  ; b ] =
    [   _1^⊤   _2^⊤; _N_1^⊤ _N_2^⊤ ][  _N_1; -_N_2 ]


and finally


    [ (_w + λ) + N_1_1_1^⊤ + N_2_2_2^⊤                                N;                             N ^⊤                                N ][  ; b ] =
    [ N_1_1 - N_2_2;     N_1 - N_2 ]

The rest of the proof follows the approach in the previous section, with_wbeing replaced by_w + λ. This proves the normal equations in eq:ridge-normal-equations-Xa correspond to regularised LDA.






§ PROOF OF LEMMA

Let _b  = λ _w be the generalised eigenvalue problem associated with a binary classification problem with unequal class means (_1≠_2) and let _w be positive definite. Then there is one non-zero eigenvalue λ = N_1 N_2/N (_1-_2)^⊤_w^-1 (_1-_2) > 0. The associated eigenvector is proportional to _w^-1(_1-_2).

Define Δ := _1 - _2 and := _w^-1Δ. Since _w is regular, the generalised eigenvalue problem can be written as an ordinary eigenvalue problem  _w^-1_b  = λ. Then using eq:Sb-simple for _b one obtains


    _w^-1_b = _w^-1(N_1 N_2/NΔΔ^⊤) = _w^-1Δ(N_1 N_2/N Δ^⊤_w^-1Δ_:=λ) = λ ,


hence  is an eigenvector of _w^-1_b with eigenvalue λ. Since _w^-1 is positive definite, λ > 0. Since _w^-1_b is of rank 1, all other eigenvalues are zero.















§ ACKNOWLEDGEMENTS


I would like to thank Richard Henson for helpful comments.










§ REFERENCES
model1-num-names