Równoległe obliczanie macierzy dużych kowariancji

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ń.$10000\times10000$$100000\times100000$

Pierwszą rzeczą jest uznanie, że możesz to zrobić za pomocą BLAS. Jeśli macierz danych jest$X = [x_1 x_2 x_3 ...] \in \mathbb{R}^{m\times n}$ (każdy $x$jest wektorem kolumny odpowiadającym jednemu pomiarowi; wiersze są próbami), możesz zapisać kowariancję jako:

$$C_{ij} = E[x_i,x_j] - E[x_i] E[x_j] = \frac{1}{n} \sum_k x_{ik} x_{jk} - \frac{1}{n^2} \left(\sum_k x_{ik} \right) \left(\sum_k x_{jk}\right)$$ Możemy to napisać jako: $$C = \frac{1}{n} X^T X - \frac{1}{n^2} (1^T X)^T (1^T X)$$ gdzie $(1^T)$ to wektor wiersza ze wszystkimi elementami 1 tak $(1^T X)$ to wektor wiersza sum kolumn $X$. Można to zapisać całkowicie jako BLAS, gdzie$X^T X$jest albo GEMM, albo jeszcze lepiej SYRK / HERK i możesz go zdobyć$(1^T X) = b$z GEMV ,$b^T b$ponownie 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 .

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);
}