Computing the Mahalanobis distance with CUDA

Accelerated Computing CUDA CUDA Programming and Performance
1

Hi,

I’m trying to modify the code that computes the Euclidian distance from: http://home.hiroshima-u.ac.jp/tamaki/study/cuda_softassign_emicp/ to compute instead the Mahalanobis distance.
Is there an example in CUDA ? even with cuBLAS implementation ?

Up to now i found the implementation of the covariance matrix calculation:

__global__ static void
updateA(int rowsA, int colsA, int pitchA,
	const float* d_Xx, const float* d_Xy, const float* d_Xz, 
	const float* d_Yx, const float* d_Yy, const float* d_Yz,
	const float* d_R, const float* d_t,
	float* d_A,
	float sigma_p2){

...

// compute covariance matrix using cublas

      // cublasSgemm (char transa, char transb, int m, int n, int k, float alpha, 
      //              const float *A, int lda, const float *B, int ldb, float beta, 
      //              float *C, int ldc)
      //   C = alpha * op(A) * op(B) + beta * C,
      //
      // m      number of rows of matrix op(A) and rows of matrix C
      // n      number of columns of matrix op(B) and number of columns of C
      // k      number of columns of matrix op(A) and number of rows of op(B) 
	  cublasSgemm('n', 't',
	              rowsA, rowsA, colsA, // int m, n (rows of A), k (cols of A) ; not op(A)
		      1.0f,         // float alpha
		      d_A, pitchA,  // const float *A, int lda
		      d_A, pitchA,  // const float *A, int lda
		      0.0f,         // float beta
                      d_Covariance, rowsA);

...

// #define Euclid(a,b,c) ((a)*(a)+(b)*(b)+(c)*(c))
//     float tmp =
//       Euclid(Xx - (R(0)*Yx + R(1)*Yy + R(2)*Yz + t(0)),
//              Xy - (R(3)*Yx + R(4)*Yy + R(5)*Yz + t(1)),
//              Xz - (R(6)*Yx + R(7)*Yy + R(8)*Yz + t(2)) );
    
//     tmp = expf(-tmp/sigma_p^2)


     float tmpX = Xx - (R(0)*Yx + R(1)*Yy + R(2)*Yz + t(0));
     float tmpY = Xy - (R(3)*Yx + R(4)*Yy + R(5)*Yz + t(1));
     float tmpZ = Xz - (R(6)*Yx + R(7)*Yy + R(8)*Yz + t(2));

    __syncthreads();

     tmpX *= tmpX;
     tmpY *= tmpY;
     tmpZ *= tmpZ;

     tmpX += tmpY;
     tmpX += tmpZ;

     tmpX /= sigma_p2;
     tmpX = expf(-tmpX);

Not sure if the covariance using cublasSgemm is correct.
How do i compute the inverse of the covariance matrix ?

Maybe there is a faster method like in OpenCV that has function call to mahalanobis distance.
Right now i can’t switch to an OpenCV implementation because of the time constraint.

Thank You,
Domi

2

In most cases you can use cuBLAS Strsm():

http://docs.nvidia.com/cuda/cublas/index.html#cublas-lt-t-gt-trsm

Also,if you are going to make a cuBLAS call from a device kernel (which I have never done before, usually call from host) , you are going to need a GPU which supports Dynamic Parallelism.

It makes more sense to call from host, does the above work correctly?

3

Thank you for the info on inverse calculation.

The code works (that’s how it is in the reference code so i did not change it), not sure if the values passed would yield the right result. I will need to cross reference with Matlab.

I’m just surprised nobody implemented the Mahalanobis distance computation using GPU, everybody is just saying how suitable this is for GPU.

Thanks,
Domi

4

@_domi would you mind to share working example? and can it be used to comparing two 3 dimensional Arrays?