Komputasi paralel dari matriks kovarians besar

9

Kita perlu menghitung matriks kovarians dengan ukuran mulai dari hingga 100000 × 100000 . Kami memiliki akses ke GPU dan cluster, kami bertanya-tanya apa pendekatan paralel terbaik untuk mempercepat perhitungan ini.10000×10000100000×100000

Buka jalannya
sumber
1
Apakah Anda mengharapkan kekhususan untuk matriks kovarian Anda? Misalnya, sejumlah besar korelasi "mendekati 0"?
Dr_Sam
tidak, kami tidak dapat mengharapkan apa pun sekarang
Buka jalan
Apa k-mu? Artinya, berapa lama masing-masing vektor data. Apakah mereka sudah nol berarti?
Max Hutchinson
tidak, mereka tidak berarti nol, mereka dapat mengambil nilai apa pun
Buka jalan
3
@flow: '' data klinis '' adalah aplikasi, tetapi bukan penggunaannya. Pertanyaan saya adalah: Misalkan Anda memiliki matriks kovarian, apa yang akan Anda lakukan dengannya (dari sudut pandang matematika)? Alasan saya bertanya adalah bahwa pada akhirnya seseorang selalu menghitung sangat sedikit darinya, dan jika ini diperhitungkan, seseorang biasanya dapat mempercepat banyak hal dengan menghindari menghitung matriks kovarian penuh sementara masih mendapatkan hasil yang diinginkan berikutnya.
Arnold Neumaier

Jawaban:

17

X=[x1x2x3...]Rm×nx

Cij=E[xi,xj]E[xi]E[xj]=1nkxikxjk1n2(kxik)(kxjk)
C=1nXTX1n2(1TX)T(1TX)
(1T)(1TX)XXTX(1TX)=bbTb

Matriks data dan hasil Anda bisa sekitar 64GB, sehingga Anda tidak akan muat pada satu node atau nilai node GPU. Untuk kluster non-GPU, Anda mungkin ingin melihat PBLAS , yang terasa seperti scalapack. Untuk GPU, perpustakaan multi-simpul belum ada di sana. Magma memiliki semacam implementasi paralel BLAS yang mendasari, tetapi mungkin tidak ramah pengguna. Saya tidak berpikir CULA melakukan multi-node, tetapi sesuatu yang harus diperhatikan . CUBLAS adalah single-node.

Saya juga menyarankan agar Anda sangat mempertimbangkan untuk mengimplementasikan paralelisme sendiri, terutama jika Anda terbiasa dengan MPI dan harus menghubungkan ini ke basis kode yang ada. Dengan begitu, Anda dapat beralih antara CPU dan GPU BLAS dengan mudah dan mulai dan akhiri dengan data di tempat yang Anda inginkan. Anda tidak perlu lebih dari beberapa panggilan MPI_ALLREDUCE .

Max Hutchinson
sumber
Terima kasih atas analisis dan daftar fungsi BLAS yang relevan. Setelah membaca jawaban Anda, saya telah menggunakan ini untuk mempercepat dan mengoptimalkan perhitungan matriks kovarians dalam versi pengembangan Scilab (www.scilab.org).
Stéphane Mottelet
E[xi,xj]E[xi]E[xj]
1

Saya menerapkan formula yang diberikan oleh @ Max Hutchinson dengan CUBlas dan Cuda Thrust dan dibandingkan dengan alat perhitungan co variance online. Tampaknya milikku menghasilkan hasil yang baik. Kode di bawah ini direncanakan untuk QDA Bayes. Jadi matriks yang diberikan mungkin mengandung lebih dari satu kelas. Jadi beberapa matriks varians dihitung. Semoga bermanfaat bagi seseorang.

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