Параллельное вычисление больших ковариационных матриц

9

Нам нужно вычислить ковариационные матрицы с размерами от до . У нас есть доступ к графическим процессорам и кластерам, мы задаемся вопросом, каков наилучший параллельный подход для ускорения этих вычислений.10000×10000100000×100000

Открой путь
источник
1
Ожидаете ли вы особенностей для вашей ковариационной матрицы? Например, большое количество «около 0» корреляций?
Dr_Sam
нет, мы не можем ожидать ничего прямо сейчас
Откройте путь
Какая у тебя к? То есть как долго каждый вектор данных. Они уже ничего не значат?
Макс Хатчинсон
нет, они не имеют нулевого значения, они могут принимать любое значение
Откройте путь
3
@flow: '' клинические данные '' - это применение, но не использование. Мой вопрос был: предположим, у вас есть ковариационная матрица, что вы собираетесь с ней делать (с математической точки зрения)? Причина, по которой я спрашиваю, состоит в том, что в итоге всегда очень мало вычисляют из этого, и, если принять это во внимание, обычно можно значительно ускорить процесс, избегая вычисления полной ковариационной матрицы, в то же время получая желаемый последующий результат.
Арнольд Ноймайер

Ответы:

17

Прежде всего, нужно признать, что вы можете сделать это с помощью BLAS. Если ваша матрица данных имеет вид (каждый - это вектор-столбец, соответствующий одному измерению; строки - это испытания), тогда вы можете написать ковариация как: Мы можем записать это как: где - вектор строки со всеми элементами 1, поэтому - вектор строки сумм столбцов . Это можно записать полностью как BLAS, гдеX=[x1x2x3...]Rm×nx

Cij=E[xi,xj]E[xi]E[xj]=1nkxikxjk1n2(kxik)(kxjk)
C=1nXTX1n2(1TX)T(1TX)
(1T)(1TX)XXTXявляется либо GEMM или, еще лучше, SYRK / Herk и вы можете получить с GEMV , с ГЭММ или SYRK / Herk снова, и prefactors с SCAL .(1TX)=bbTb

Ваши матрицы данных и результатов могут составлять около 64 ГБ, поэтому вы не собираетесь помещать их на один узел или на графические процессоры. Для кластера без графического процессора вы можете посмотреть на PBLAS , который выглядит как скальпак. Для графических процессоров многоузловые библиотеки еще не совсем там. Magma имеет своего рода параллельную реализацию BLAS, но она не может быть удобной для пользователя. Я не думаю, что CULA еще работает с несколькими узлами, но это то, за чем нужно следить. CUBLAS является одноузловым .

Я бы также посоветовал вам настоятельно рассмотреть возможность реализации параллелизма самостоятельно, особенно если вы знакомы с MPI и должны подключить это к существующей кодовой базе. Таким образом, вы можете легко переключаться между процессором и GPU BLAS и начинать и заканчивать данными именно там, где вам нужно. Вам не нужно больше, чем несколько вызовов MPI_ALLREDUCE .

Макс Хатчинсон
источник
Спасибо за ваш анализ и список соответствующих функций BLAS. После прочтения вашего ответа я использовал их для ускорения и оптимизации вычисления ковариационной матрицы в версии Scilab для разработчиков (www.scilab.org).
Стефан Моттлет
Однако следует помнить, что использование этого способа вычисления ковариации может привести к катастрофической отмене, когда близко к , см., Например, en.wikipedia.org/wiki/…E[xi,xj]E[xi]E[xj]
Stéphane Mottelet
1

Я реализовал формулу, заданную @Max Hutchinson с CUBlas и Cuda Thrust, и сравнил ее с онлайн-инструментами вычисления дисперсии. Кажется, у меня хорошие результаты. Код ниже планируется QDA Байеса. Таким образом, данная матрица может содержать более одного класса. Таким образом, вычисляется несколько матриц. Надеюсь это кому-нибудь пригодится.

//! 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);
}
Кадир Эрдем Демир
источник