Cálculo paralelo de grandes matrizes de covariância

9

Precisamos calcular matrizes de covariância com tamanhos que variam de a . Temos acesso a GPUs e clusters, nos perguntamos qual é a melhor abordagem paralela para acelerar esses cálculos.100000 × 10000010000×10000100000×100000

Abra o caminho
fonte
11
Você espera particularidades para sua matriz de covariância? Por exemplo, grande número de correlações "próximas de 0"?
21813 Dr_Sam # 03
não, não podemos esperar nada agora
Abra o caminho
Qual é o seu k? Ou seja, quanto tempo dura cada vetor de dados. Eles já têm média zero?
Max Hutchinson
Não, eles não são de média zero, eles podem tomar qualquer valor
Abrir o caminho
3
@flow: '' dados clínicos '' é a aplicação, mas não o uso. Minha pergunta era: suponha que você tenha a matriz de covariância, o que você fará com ela (do ponto de vista matemático)? A razão pela qual pergunto é que, no final, sempre se calcula muito pouco e, se isso é levado em consideração, geralmente é possível acelerar tremendamente as coisas, evitando calcular a matriz de covariância completa e ainda assim obter o resultado subsequente desejado.
Arnold Neumaier 11/03/2013

Respostas:

17

A primeira coisa é reconhecer que você pode fazer isso usando o BLAS. Se sua matriz de dados é (cada é um vetor de coluna correspondente a uma medida; as linhas são tentativas), você pode escrever a covariância como: Podemos escrever isso como: onde é o vetor de linha com todos os elementos 1, de modo que é um vetor de linha da soma da coluna de . Isso pode ser escrito completamente como BLAS, onde x C i j = E [ x i , x j ] - E [ x i ] E [ x j ] = 1X=[x1 1x2x3...]Rm×nx

CEuj=E[xEu,xj]-E[xEu]E[xj]=1 1nkxEukxjk-1 1n2(kxEuk)(kxjk)
C=1 1nXTX-1 1n2(1 1TX)T(1 1TX)
(1 1T)(1 1TX)XXTXé um GEMM ou, melhor ainda, um SYRK / HERK e você pode obter com um GEMV , Tb com GEMM ou SYRK / HERK novamente e os prefatores com SCAL .(1 1TX)=bbTb

Suas matrizes de dados e resultados podem ter cerca de 64 GB, portanto, você não vai caber em um único nó ou nas GPUs de um nó. Para um cluster sem GPU, convém consultar o PBLAS , que parece um escalapack. Para GPUs, as bibliotecas de vários nós ainda não estão lá. O Magma possui algum tipo de implementação paralela do BLAS, mas pode não ser fácil de usar. Eu não acho que o CULA faça vários nós ainda, mas é algo para ficar de olho. CUBLAS é um nó.

Eu também sugiro que você considere fortemente implementar o paralelismo, especialmente se você estiver familiarizado com o MPI e precisar conectá-lo a uma base de código existente. Dessa forma, você pode alternar entre CPU e GPU BLAS facilmente e começar e terminar com os dados exatamente onde deseja. Você não precisa de mais do que algumas chamadas MPI_ALLREDUCE .

Max Hutchinson
fonte
Obrigado por sua análise e lista de funções BLAS relevantes. Depois de ler sua resposta, eu as usei para acelerar e otimizar o cálculo da matriz de covariância na versão de desenvolvimento do Scilab (www.scilab.org).
Stéphane Mottelet 8/04/19
No entanto, esteja avisado de que, usando essa maneira de calcular, a covariância está sujeita a cancelamento catastrófico quando estiver próximo de , consulte, por exemplo, en.wikipedia.org/wiki/…E[xEu,xj]E[xEu]E[xj]
Stéphane Mottelet
1

Eu implementei a fórmula dada por @Max Hutchinson com CUBlas e Cuda Thrust e comparei com as ferramentas de cálculo de variação on-line. Parece o meu produzindo bons resultados. O código abaixo planejado para QDA Bayes. Portanto, a matriz fornecida pode conter mais de uma classe. Portanto, várias matrizes de co-variação são calculadas. Espero que seja útil para alguém.

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