Parallele Berechnung großer Kovarianzmatrizen

9

Wir müssen Kovarianzmatrizen mit Größen von bis 100000 × 100000 berechnen . Wir haben Zugriff auf GPUs und Cluster und fragen uns, was der beste parallele Ansatz ist, um diese Berechnungen zu beschleunigen.10000×10000100000×100000

Mach den Weg frei
quelle
1
Erwarten Sie Besonderheiten für Ihre Kovarianzmatrix? Zum Beispiel eine große Anzahl von "nahe 0" -Korrelationen?
Dr_Sam
Nein, wir können im Moment nichts erwarten
Öffnen Sie den Weg
Was ist dein k? Das heißt, wie lang ist jeder Datenvektor. Sind sie schon Null-Mittelwert?
Max Hutchinson
Nein, sie sind nicht Null-Mittelwert, sie können jeden Wert annehmen
Öffnen Sie den Weg
3
@flow: '' klinische Daten '' ist die Anwendung, aber nicht die Verwendung. Meine Frage war: Angenommen, Sie haben die Kovarianzmatrix, was werden Sie damit machen (aus mathematischer Sicht)? Der Grund, den ich frage, ist, dass man am Ende immer sehr wenig daraus berechnet, und wenn dies berücksichtigt wird, kann man die Dinge normalerweise enorm beschleunigen, indem man vermeidet, die vollständige Kovarianzmatrix zu berechnen, während man immer noch das gewünschte nachfolgende Ergebnis erhält.
Arnold Neumaier

Antworten:

17

Als erstes müssen Sie erkennen, dass Sie dies mit BLAS tun können. Wenn Ihre Datenmatrix (jedes x ist ein Spaltenvektor, der einer Messung entspricht; Zeilen sind Versuche), dann können Sie die Kovarianz wie folgt schreiben: C i j = E [ x i , x j ] - E [ x i ] E. [ x j ] = 1X=[x1x2x3...]Rm×nx Wir können dies schreiben als: C=1

Cij=E[xi,xj]E[xi]E[xj]=1nkxikxjk1n2(kxik)(kxjk)
wobei(1T)der Zeilenvektor mit allen Elementen 1 ist, so dass(1TX)ein Zeilenvektor der Spaltensummen vonX ist. Dies kann vollständig als BLAS geschrieben werden, wobeiXTXentweder einGEMModer besser noch einSYRK/HERK istund Siemit einemGEMV(1TX)=berhaltenkönnen.
C=1nXTX1n2(1TX)T(1TX)
(1T)(1TX)XXTX(1TX)=b mit GEMM oder SYRK / HERK und die Vorfaktoren mitSCAL.bTb

Ihre Daten- und Ergebnismatrizen können etwa 64 GB groß sein, sodass Sie nicht auf einen einzelnen Knoten oder die GPUs eines Knotens passen. Für einen Nicht-GPU-Cluster sollten Sie sich PBLAS ansehen , das sich wie ein Scalapack anfühlt. Für GPUs sind Bibliotheken mit mehreren Knoten noch nicht ganz vorhanden. Magma hat eine zugrunde liegende parallele BLAS-Implementierung, die jedoch möglicherweise nicht benutzerfreundlich ist. Ich denke, CULA macht noch keine Multi-Node, aber es ist etwas, das man im Auge behalten sollte. CUBLAS ist ein Einzelknoten .

Ich würde auch vorschlagen, dass Sie die Parallelität unbedingt selbst implementieren, insbesondere wenn Sie mit MPI vertraut sind und dies in eine vorhandene Codebasis einbinden müssen. Auf diese Weise können Sie einfach zwischen CPU und GPU BLAS wechseln und mit den Daten genau dort beginnen und enden, wo Sie sie möchten. Sie sollten nicht mehr als ein paar MPI_ALLREDUCE- Aufrufe benötigen .

Max Hutchinson
quelle
Vielen Dank für Ihre Analyse und Liste der relevanten BLAS-Funktionen. Nachdem ich Ihre Antwort gelesen habe, habe ich diese verwendet, um die Berechnung der Kovarianzmatrix in der Entwicklungsversion von Scilab (www.scilab.org) zu beschleunigen und zu optimieren.
Stéphane Mottelet
E[xi,xj]E[xi]E[xj]
1

Ich habe die von @Max Hutchinson gegebene Formel mit CUBlas und Cuda Thrust implementiert und mit Online-Tools zur Berechnung der Varianz verglichen. Es scheint, als würde ich gute Ergebnisse erzielen. Der folgende Code ist für QDA Bayes geplant. Die angegebene Matrix kann also mehr als eine Klasse enthalten. Es werden also mehrere Co-Varianz-Matrizen berechnet. Ich hoffe, es wird für jemanden nützlich sein.

//! 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
quelle