Calcul parallèle de grandes matrices de covariance

9

Nous devons calculer des matrices de covariance avec des tailles allant de à . Nous avons accès aux GPU et aux clusters, nous nous demandons quelle est la meilleure approche parallèle pour accélérer ces calculs.100000 × 10000010000×10000100000×100000

Ouvrez la voie
la source
1
Vous attendez-vous à des particularités pour votre matrice de covariance? Par exemple, un grand nombre de corrélations "proches de 0"?
Dr_Sam
non, nous ne pouvons rien attendre pour l'instant
Ouvrez la voie
Quel est ton k? Autrement dit, quelle est la longueur de chaque vecteur de données. Sont-ils déjà à moyenne nulle?
Max Hutchinson
non, ils ne sont pas à moyenne nulle, ils peuvent prendre n'importe quelle valeur
Ouvrez la voie
3
@flow: '' données cliniques '' est l'application, mais pas l'utilisation. Ma question était: supposons que vous ayez la matrice de covariance, qu'allez-vous en faire (d'un point de vue mathématique)? La raison pour laquelle je demande, c'est qu'en fin de compte, on en calcule toujours très peu, et si cela est pris en compte, on peut généralement accélérer considérablement les choses en évitant de calculer la matrice de covariance complète tout en obtenant le résultat souhaité.
Arnold Neumaier

Réponses:

17

La première chose est de reconnaître que vous pouvez le faire en utilisant BLAS. Si votre matrice de données est (chaque est un vecteur de colonne correspondant à une mesure; les lignes sont des essais), alors vous pouvez écrire la covariance comme: Nous pouvons écrire ceci comme: où est le vecteur ligne avec tous les éléments 1 donc est un vecteur ligne de la somme des colonnes de . Cela peut être écrit complètement comme BLAS, où x C i j = E [ x i , x j ] - E [ x i ] E [ x j ] = 1X=[X1X2X3...]Rm×nX

Cjej=E[Xje,Xj]-E[Xje]E[Xj]=1nkXjekXjk-1n2(kXjek)(kXjk)
(1T)(1TX)XXTX(1TX)=bbTb
C=1nXTX-1n2(1TX)T(1TX)
(1T)(1TX)XXTXest soit un GEMM ou, mieux encore, un SYRK / HERK et vous pouvez obtenir le avec un GEMV , avec GEMM ou SYRK / HERK à nouveau, et les préfacteurs avec SCAL .(1TX)=bbTb

Vos données et matrices de résultats peuvent faire environ 64 Go, donc vous n'allez pas tenir sur un seul nœud ou sur la valeur d'un nœud de GPU. Pour un cluster non GPU, vous voudrez peut-être regarder PBLAS , qui ressemble à un scalapack. Pour les GPU, les bibliothèques multi-nœuds n'en sont pas encore là. Magma a une sorte d'implémentation BLAS parallèle sous-jacente, mais elle peut ne pas être conviviale. Je ne pense pas que CULA fasse encore plusieurs nœuds, mais c'est quelque chose à surveiller. CUBLAS est à nœud unique.

Je suggérerais également que vous envisagiez fortement d'implémenter le parallélisme vous-même, surtout si vous êtes familier avec MPI et devez le connecter à une base de code existante. De cette façon, vous pouvez basculer facilement entre CPU et GPU BLAS et commencer et terminer avec les données exactement où vous le souhaitez. Vous ne devriez pas avoir besoin de plus de quelques appels MPI_ALLREDUCE .

Max Hutchinson
la source
Merci pour votre analyse et la liste des fonctions BLAS pertinentes. Après avoir lu votre réponse, je les ai utilisées pour accélérer et optimiser le calcul de la matrice de covariance dans la version de développement de Scilab (www.scilab.org).
Stéphane Mottelet
E[Xje,Xj]E[Xje]E[Xj]
1

J'ai implémenté la formule donnée par @Max Hutchinson avec CUBlas et Cuda Thrust et comparé avec des outils de calcul de co-variance en ligne. Il semble que le mien produise de bons résultats. Le code ci-dessous prévu pour QDA Bayes. La matrice donnée peut donc contenir plus d'une classe. Ainsi, plusieurs matrices de co-variance sont calculées. J'espère que ce sera utile pour quelqu'un.

//! Calculates one or more than one coVarianceMatrix given data.
//  There can be many classes since many covariance matrixes.
/*!
    \param inMatrix This vector contains matrix data in major storage. 
    Forexample if inMatrix=[1 2 3 4 5 6] and trialSizes=[2] this means matrix we will work on a matrix like :
        |1 4 |
        |2 5 |
        |3 6 | -> 2 Trials, 3 Features. Columns contains feature rows contains trials (samples)
    \param trialSizes There can be many classes since many covariance matrixes. Samples from all classes will be given with inMatrix.
    But we need to know how many trials(samples) we have for each class. 
    For example if inMatrix=[1 2 3 4 5 6 7 8 9 10 11 12] and trialSizes=[2,2] 
    this means matrix we will work on a matrix like :
        |1 4 |  |7 10 |
        |2 5 |  |8 11 |
        |3 6 |  |9 12 |  --> Total number of trials(samples which is total rowCount) 2 + 2 = 4 , 
                             So colSize = inMatrix.size()/4 = 3(feature vector size)
                         --> There is two element in trialSize vec so each vector has to samples
*/
void multiQDACovianceCalculator(std::vector<float>& inMatrix, std::vector<int>& trialSizes)
{
    cublasHandle_t handle; // CUBLAS context
    int classCount = trialSizes.size();
    int rowSize = std::accumulate(trialSizes.begin(), trialSizes.end(), 0);
    int dimensionSize = inMatrix.size() / rowSize;
    float alpha = 1.0f;
    float beta = 0.0f; // bet =1

    thrust::device_vector<float> d_cov1(dimensionSize * dimensionSize);
    thrust::device_vector<float> d_cov2(dimensionSize * dimensionSize);
    thrust::device_vector<float> d_covResult(dimensionSize * dimensionSize);

    thrust::device_vector<float> d_wholeMatrix(inMatrix);
    thrust::device_vector<float> d_meansVec(dimensionSize); // rowVec of means of trials
    float *meanVecPtr = thrust::raw_pointer_cast(d_meansVec.data());
    float *device2DMatrixPtr = thrust::raw_pointer_cast(d_wholeMatrix.data());
    auto maxTrialNumber = *std::max_element(trialSizes.begin(), trialSizes.end());
    thrust::device_vector<float> deviceVector(maxTrialNumber, 1.0f);

    cublasCreate(&handle);
    // Inside of for loop  one covariance matrix calculated each time
    for (int i = 0; i < trialSizes.size(); i++)
    {
        // X*transpose(X) / N
        alpha = 1.0f / trialSizes[i];
        cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, dimensionSize, dimensionSize, trialSizes[i], &alpha,
            device2DMatrixPtr, dimensionSize, device2DMatrixPtr, dimensionSize, &beta,
            thrust::raw_pointer_cast(d_cov1.data()), dimensionSize);

        // Mean vector of each column
        alpha = 1.0f;
        cublasSgemv(handle, CUBLAS_OP_N, dimensionSize, trialSizes[i], &alpha, device2DMatrixPtr,
            dimensionSize, thrust::raw_pointer_cast(deviceVector.data()), 1, &beta, meanVecPtr, 1);

        // MeanVec * transpose(MeanVec) / N*N
        alpha = 1.0f / (trialSizes[i] * trialSizes[i]);
        cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, dimensionSize, dimensionSize, 1, &alpha,
            meanVecPtr, 1, meanVecPtr, 1, &beta,
            thrust::raw_pointer_cast(d_cov2.data()), dimensionSize);

        alpha = 1.0f;
        beta = -1.0f;
        //  (X*transpose(X) / N) -  (MeanVec * transpose(MeanVec) / N*N)
        cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, dimensionSize, dimensionSize, &alpha,
            thrust::raw_pointer_cast(d_cov1.data()), dimensionSize, &beta, thrust::raw_pointer_cast(d_cov2.data()), 
            dimensionSize, thrust::raw_pointer_cast(d_covResult.data()), dimensionSize);

        // Go to other class and calculate its covarianceMatrix
        device2DMatrixPtr += trialSizes[i] * dimensionSize;
    }
    printVector(d_covResult);
    cublasDestroy(handle);
}
Kadir Erdem Demir
la source