How to compute pairwise Euclidean distances using IPP?

[minne]

I want to compute the pairwise squared Euclidean distances between all the rows of an NxD matrix and store these distances in an NxN matrix. In plain C++ code, that would look something like this:

for(int i = 0; i smaller than n; i++) {

for(int j = i + 1; j smaller than n; j++) {

double val = 0.0;

for(int k = 0; k
val += pow((float) *(data + (i * (data_step / sizeof(Ipp32f))) + k) - (float) *(data + (j * (data_step / sizeof(Ipp32f))) + k), (float) 2.0);

}

*(P + (i * (P_step / sizeof(Ipp32f))) + j) = *(P + (j * (P_step / sizeof(Ipp32f))) + i) = (Ipp32f) val;

}

}

Herein, the data is my NxD matrix with stride data_step, and P is the NxN matrix in which I store the squared Euclidean distances.

I've been trying to think about a smart way to do this with IPP (e.g., by using d(A,B) = sum(A.^2) + sum(B.^2) - 2*sum(A.*B)), but they all involve the implementation of loops. Does anyone have a suggestion how I can do this computation efficiently using IPP?

[Intel_C_Intel]

Hi,

It is possible to apply one of ippsLogGauss_IdVar_* function (ippsr library). The formula slightly differs (for val==0 you get the negative Euclidean square distance.

32f, 16s to 32s and 16s to 32f data types are supported.

Data step in these functions is in data elements, not in bytes

Alexander

[minne]

Thanks, but that only allows me to compute one row of the pairwise distance matrix at a time right? I've tried that, but it doesn't really make things any faster (since you don't have the advantage of processing a large chunk of data at a time)...

[minne]

The best I came up with until now is the following (assuming you have an NxD dataset called data):

=========================
Ipp64f* dataSqr = ippsMalloc_64f (n * d);
Ipp64f* dataSqrSums = ippsMalloc_64f (n);
ippsSqr_64f (data, dataSqr, n * d);
ippsSumRow_64f_D2 (dataSqr, d, d, dataSqrSums, n);
ippmMul_mt_64f (data, d * sizeof(Ipp64f), sizeof(Ipp64f), d, n, data, d * sizeof(Ipp64f), sizeof(Ipp64f), d, n, P, n * sizeof(Ipp64f), sizeof(Ipp64f));
ippsMulC_64f_I ((Ipp64f) -2.0, P, n * n);
for(int i = 0; i smaller than n; i++) {
for(int j = 0; j smaller than n; j++) {
*(P + i * n + j) += *(dataSqrSums + i) + *(dataSqrSums + j);
}
}
ippsFree(dataSqr);
ippsFree(dataSqrSums);
==========================

However, it has three disadvantages:
* It uses ippm, which I think is mainly optimized for small matrices (in my case, N=5000 and D=[1,10])
* I still need a loop to add the L2-norms of the data vectors
* I don't exploit the fact that the pairwise distance matrix is symmetric (I compute everything twice)

Any ideas anyone?