Równoległe obliczanie macierzy dużych kowariancji

9

Musimy obliczyć macierze kowariancji o rozmiarach od do . Mamy dostęp do GPU i klastrów, zastanawiamy się, jakie jest najlepsze równoległe podejście do przyspieszenia tych obliczeń.dziesięć tysięcy×dziesięć tysięcy100000×100000

Otwórz drogę
źródło
1
Czy oczekujesz szczególnych cech macierzy kowariancji? Na przykład duża liczba korelacji „w pobliżu 0”?
Dr_Sam
nie, nie możemy się teraz niczego spodziewać
otwórz drogę
Jaki jest twój k To znaczy, jak długi jest każdy wektor danych. Czy są już zerowe?
Max Hutchinson
nie, nie są one równe zeru, mogą przyjąć dowolną wartość
otwórz sposób
3
@flow: „dane kliniczne” to appl8ication, ale nie zastosowanie. Moje pytanie brzmiało: Załóżmy, że masz macierz kowariancji, co zamierzasz z nią zrobić (z matematycznego punktu widzenia)? Pytam dlatego, że ostatecznie zawsze obliczamy z tego bardzo mało, a jeśli się to weźmie pod uwagę, zwykle można ogromnie przyspieszyć, unikając obliczenia pełnej macierzy kowariancji, a jednocześnie uzyskać pożądany wynik.
Arnold Neumaier

Odpowiedzi:

17

Pierwszą rzeczą jest uznanie, że możesz to zrobić za pomocą BLAS. Jeśli macierz danych jestX=[x1x2x3...]Rm×n (każdy xjest wektorem kolumny odpowiadającym jednemu pomiarowi; wiersze są próbami), możesz zapisać kowariancję jako:

dojajot=mi[xja,xjot]-mi[xja]mi[xjot]=1nkxjakxjotk-1n2)(kxjak)(kxjotk)
Możemy to napisać jako:
do=1nXT.X-1n2)(1T.X)T.(1T.X)
gdzie (1T.) to wektor wiersza ze wszystkimi elementami 1 tak (1T.X) to wektor wiersza sum kolumn X. Można to zapisać całkowicie jako BLAS, gdzieXT.Xjest albo GEMM, albo jeszcze lepiej SYRK / HERK i możesz go zdobyć(1T.X)=bz GEMV ,bT.bponownie z GEMM lub SYRK / HERK, a prefaktory z SCAL .

Twoje dane i macierze wyników mogą mieć około 64 GB, więc nie zmieścisz się na pojedynczym węźle ani na GPU o wartości jednego węzła. W przypadku klastra bez GPU warto przyjrzeć się PBLAS , który wydaje się być scalapack. W przypadku układów GPU biblioteki wielu węzłów jeszcze nie są dostępne. Magma ma jakąś podstawową równoległą implementację BLAS, ale może nie być przyjazna dla użytkownika. Nie sądzę, że CULA ma jeszcze wiele węzłów, ale warto mieć na to oko. CUBLAS jest jednowęzłowy.

Sugerowałbym również, abyś mocno zastanowił się nad implementacją równoległości, szczególnie jeśli znasz MPI i musisz podłączyć go do istniejącej bazy kodu. W ten sposób możesz łatwo przełączać się między procesorem a GPU BLAS i zaczynać i kończyć danymi dokładnie tam, gdzie chcesz. Nie powinieneś potrzebować więcej niż kilku wywołań MPI_ALLREDUCE .

Max Hutchinson
źródło
Dziękujemy za analizę i listę odpowiednich funkcji BLAS. Po przeczytaniu twojej odpowiedzi wykorzystałem je do przyspieszenia i optymalizacji obliczeń macierzy kowariancji w rozwojowej wersji Scilab (www.scilab.org).
Stéphane Mottelet,
Należy jednak pamiętać, że stosowanie tego sposobu obliczania kowariancja podlega katastrofalnej anulowaniu, gdy mi[xja,xjot] jest blisko do mi[xja]mi[xjot], patrz np. en.wikipedia.org/wiki/…
Stéphane Mottelet
1

Zaimplementowałem formułę podaną przez @Max Hutchinson z CUBlas i Cuda Thrust i porównałem z narzędziami do obliczania wariancji online. Wydaje się, że mój przynosi dobre wyniki. Poniższy kod planowany jest dla QDA Bayes. Podana macierz może zawierać więcej niż jedną klasę. Tak więc oblicza się wiele macierzy wariancji co. Mam nadzieję, że przyda się komuś.

//! 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
źródło