推力矢量距離計算

Question

考慮以下數據集和質心。有8個維度的7個人和兩個手段。他們存儲行主要訂單。

short dim = 8; 
float centroids[] = { 
    0.223, 0.002, 0.223, 0.412, 0.334, 0.532, 0.244, 0.612, 
    0.742, 0.812, 0.817, 0.353, 0.325, 0.452, 0.837, 0.441 
}; 
float data[] = { 
    0.314, 0.504, 0.030, 0.215, 0.647, 0.045, 0.443, 0.325, 
    0.731, 0.354, 0.696, 0.604, 0.954, 0.673, 0.625, 0.744, 
    0.615, 0.936, 0.045, 0.779, 0.169, 0.589, 0.303, 0.869, 
    0.275, 0.406, 0.003, 0.763, 0.471, 0.748, 0.230, 0.769, 
    0.903, 0.489, 0.135, 0.599, 0.094, 0.088, 0.272, 0.719, 
    0.112, 0.448, 0.809, 0.157, 0.227, 0.978, 0.747, 0.530, 
    0.908, 0.121, 0.321, 0.911, 0.884, 0.792, 0.658, 0.114 
}; 

我想計算每個歐幾里德距離。 C1 - D1，C1 - D2 .... CPU上，我會做：

float dist = 0.0, dist_sqrt; 
for(int i = 0; i < 2; i++) 
    for(int j = 0; j < 7; j++) 
    { 
     float dist_sum = 0.0; 
     for(int k = 0; k < dim; k++) 
     { 
      dist = centroids[i * dim + k] - data[j * dim + k]; 
      dist_sum += dist * dist; 
     } 
     dist_sqrt = sqrt(dist_sum); 
     // do something with the distance 
     std::cout << dist_sqrt << std::endl; 

    } 

有沒有建在推力矢量距離計算的解決方案？

Answer 1

它可以在推力。解釋將如何涉及，代碼密集。

從關鍵的觀察開始，核心操作可以通過轉換的減少來完成。推力變換操作用於執行每個結果的矢量（單個質心）和平方的元素減法，並且減法將結果相加在一起以產生歐幾里得距離的平方。該操作的起點是thrust::reduce_by_key，但它涉及到正確地將數據呈現給reduce_by_key。

最終的結果是通過從上面得到每個結果的平方根產生的，我們可以使用普通的thrust::transform。

以上是對所有工作進行的唯一兩行推力代碼的總結說明。但是，第一行相當複雜。爲了利用並行性，我採用的方法是實際上按順序「佈局」必要的向量，以呈現給reduce_by_key。舉一個簡單的例子，假設我們有2個質心和4個人，假設我們的尺寸爲2

centroid 0: C00 C01 
centroid 1: C10 C11 
individ 0: I00 I01 
individ 1: I10 I11 
individ 2: I20 I21 
individ 3: I30 I31 

我們可以在「佈局」的載體是這樣的：

C00 C01 C00 C01 C00 C01 C00 C01 C10 C11 C10 C11 C10 C11 C10 C11 
I00 I01 I10 I11 I20 I21 I30 I31 I00 I01 I10 I11 I20 I21 I30 I31 

爲了方便reduce_by_key，我們還需要將生成鍵值來描繪載體：

0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 

以上數據「下崗」數據集是非常大的，我們不希望引發一個儲存d檢索成本，因此我們將使用推力的收集fancy iterators來生成這些「即時」。這是事情變得非常密集的地方。考慮到上述策略，我們將使用thrust::reduce_by_key來完成這項工作。我們將創建一個提供給transform_iterator的自定義函數，用於對I和C向量進行減法（和平方），這些向量將被壓縮在一起以達到此目的。矢量的「佈局」將使用具有附加定製索引創建函數的置換迭代器實時創建，以幫助使用I和C中的每一箇中的複製模式。

因此，從「內向外」工作，步驟的順序如下：

1. 兩個I（data）和C（centr）使用counting_iterator內部具有自定義索引函子組合的transform_iterator來產生我們需要的索引序列。

2. 使用在步驟1和基座I和C矢量創建的索引序列，經由permutation_iterator（每個佈局的矢量）虛擬地「鋪陳」的載體。

3. 將2個「佈局」虛擬I和C向量拉到一起，創建一個<float, float>元組向量（虛擬）。

4. 採取從步驟3 zip_iterator，並將其與以transform_iterator

5. 使用另一transform_iterator定製距離計算函子（(I-C)^2）結合，與自定義密鑰產生算符組合counting_iterator，以產生鍵序列（虛擬）

6. 將步驟4和5中的迭代器傳遞給reduce_by_key作爲待減少的輸入（鍵，值）。 reduce_by_key的輸出向量也是鍵和值。我們不需要鑰匙，所以我們將使用discard_iterator來轉儲這些鑰匙。我們將保存的值。

上述步驟都是用單行推力碼完成的。

下面是示出上述的碼：

#include <iostream> 
#include <thrust/device_vector.h> 
#include <thrust/host_vector.h> 
#include <thrust/reduce.h> 
#include <thrust/iterator/transform_iterator.h> 
#include <thrust/iterator/counting_iterator.h> 
#include <thrust/iterator/permutation_iterator.h> 
#include <thrust/iterator/zip_iterator.h> 
#include <thrust/iterator/discard_iterator.h> 
#include <thrust/copy.h> 
#include <math.h> 

#include <time.h> 
#include <sys/time.h> 
#include <stdlib.h> 

#define MAX_DATA 100000000 
#define MAX_CENT 5000 
#define TOL 0.001 

unsigned long long dtime_usec(unsigned long long prev){ 
#define USECPSEC 1000000ULL 
    timeval tv1; 
    gettimeofday(&tv1,0); 
    return ((tv1.tv_sec * USECPSEC)+tv1.tv_usec) - prev; 
} 

unsigned verify(float *d1, float *d2, int len){ 
    unsigned pass = 1; 
    for (int i = 0; i < len; i++) 
    if (fabsf(d1[i] - d2[i]) > TOL){ 
     std::cout << "mismatch at: " << i << " val1: " << d1[i] << " val2: " << d2[i] << std::endl; 
     pass = 0; 
     break;} 
    return pass; 
} 
void eucl_dist_cpu(const float *centroids, const float *data, float *rdist, int num_centroids, int dim, int num_data, int print){ 

    int out_idx = 0; 
    float dist, dist_sqrt; 
    for(int i = 0; i < num_centroids; i++) 
    for(int j = 0; j < num_data; j++) 
    { 
     float dist_sum = 0.0; 
     for(int k = 0; k < dim; k++) 
     { 
      dist = centroids[i * dim + k] - data[j * dim + k]; 
      dist_sum += dist * dist; 
     } 
     dist_sqrt = sqrt(dist_sum); 
     // do something with the distance 
     rdist[out_idx++] = dist_sqrt; 
     if (print) std::cout << dist_sqrt << ", "; 

    } 
    if (print) std::cout << std::endl; 
} 


struct dkeygen : public thrust::unary_function<int, int> 
{ 
    int dim; 
    int numd; 

    dkeygen(const int _dim, const int _numd) : dim(_dim), numd(_numd) {}; 

    __host__ __device__ int operator()(const int val) const { 
    return (val/dim); 
    } 
}; 

typedef thrust::tuple<float, float> mytuple; 
struct my_dist : public thrust::unary_function<mytuple, float> 
{ 
    __host__ __device__ float operator()(const mytuple &my_tuple) const { 
    float temp = thrust::get<0>(my_tuple) - thrust::get<1>(my_tuple); 
    return temp*temp; 
    } 
}; 


struct d_idx : public thrust::unary_function<int, int> 
{ 
    int dim; 
    int numd; 

    d_idx(int _dim, int _numd) : dim(_dim), numd(_numd) {}; 

    __host__ __device__ int operator()(const int val) const { 
    return (val % (dim*numd)); 
    } 
}; 

struct c_idx : public thrust::unary_function<int, int> 
{ 
    int dim; 
    int numd; 

    c_idx(int _dim, int _numd) : dim(_dim), numd(_numd) {}; 

    __host__ __device__ int operator()(const int val) const { 
    return (val % dim) + (dim * (val/(dim*numd))); 
    } 
}; 

struct my_sqrt : public thrust::unary_function<float, float> 
{ 
    __host__ __device__ float operator()(const float val) const { 
    return sqrtf(val); 
    } 
}; 


unsigned long long eucl_dist_thrust(thrust::host_vector<float> &centroids, thrust::host_vector<float> &data, thrust::host_vector<float> &dist, int num_centroids, int dim, int num_data, int print){ 

    thrust::device_vector<float> d_data = data; 
    thrust::device_vector<float> d_centr = centroids; 
    thrust::device_vector<float> values_out(num_centroids*num_data); 

    unsigned long long compute_time = dtime_usec(0); 
    thrust::reduce_by_key(thrust::make_transform_iterator(thrust::make_counting_iterator<int>(0), dkeygen(dim, num_data)), thrust::make_transform_iterator(thrust::make_counting_iterator<int>(dim*num_data*num_centroids), dkeygen(dim, num_data)),thrust::make_transform_iterator(thrust::make_zip_iterator(thrust::make_tuple(thrust::make_permutation_iterator(d_centr.begin(), thrust::make_transform_iterator(thrust::make_counting_iterator<int>(0), c_idx(dim, num_data))), thrust::make_permutation_iterator(d_data.begin(), thrust::make_transform_iterator(thrust::make_counting_iterator<int>(0), d_idx(dim, num_data))))), my_dist()), thrust::make_discard_iterator(), values_out.begin()); 
    thrust::transform(values_out.begin(), values_out.end(), values_out.begin(), my_sqrt()); 
    cudaDeviceSynchronize(); 
    compute_time = dtime_usec(compute_time); 

    if (print){ 
    thrust::copy(values_out.begin(), values_out.end(), std::ostream_iterator<float>(std::cout, ", ")); 
    std::cout << std::endl; 
    } 
    thrust::copy(values_out.begin(), values_out.end(), dist.begin()); 
    return compute_time; 
} 


int main(int argc, char *argv[]){ 

    int dim = 8; 
    int num_centroids = 2; 
    float centroids[] = { 
    0.223, 0.002, 0.223, 0.412, 0.334, 0.532, 0.244, 0.612, 
    0.742, 0.812, 0.817, 0.353, 0.325, 0.452, 0.837, 0.441 
    }; 
    int num_data = 8; 
    float data[] = { 
    0.314, 0.504, 0.030, 0.215, 0.647, 0.045, 0.443, 0.325, 
    0.731, 0.354, 0.696, 0.604, 0.954, 0.673, 0.625, 0.744, 
    0.615, 0.936, 0.045, 0.779, 0.169, 0.589, 0.303, 0.869, 
    0.275, 0.406, 0.003, 0.763, 0.471, 0.748, 0.230, 0.769, 
    0.903, 0.489, 0.135, 0.599, 0.094, 0.088, 0.272, 0.719, 
    0.112, 0.448, 0.809, 0.157, 0.227, 0.978, 0.747, 0.530, 
    0.908, 0.121, 0.321, 0.911, 0.884, 0.792, 0.658, 0.114, 
    0.721, 0.555, 0.979, 0.412, 0.007, 0.501, 0.844, 0.234 
    }; 
    std::cout << "cpu results: " << std::endl; 
    float dist[num_data*num_centroids]; 
    eucl_dist_cpu(centroids, data, dist, num_centroids, dim, num_data, 1); 

    thrust::host_vector<float> h_data(data, data + (sizeof(data)/sizeof(float))); 
    thrust::host_vector<float> h_centr(centroids, centroids + (sizeof(centroids)/sizeof(float))); 
    thrust::host_vector<float> h_dist(num_centroids*num_data); 

    std::cout << "gpu results: " << std::endl; 
    eucl_dist_thrust(h_centr, h_data, h_dist, num_centroids, dim, num_data, 1); 

    float *data2, *centroids2, *dist2; 
    num_centroids = 10; 
    num_data = 1000000; 

    if (argc > 2) { 
    num_centroids = atoi(argv[1]); 
    num_data = atoi(argv[2]); 
    if ((num_centroids < 1) || (num_centroids > MAX_CENT)) {std::cout << "Num centroids out of range" << std::endl; return 1;} 
    if ((num_data < 1) || (num_data > MAX_DATA)) {std::cout << "Num data out of range" << std::endl; return 1;} 
    if (num_data * dim * num_centroids > 2000000000) {std::cout << "data set out of range" << std::endl; return 1;}} 
    std::cout << "Num Data: " << num_data << std::endl; 
    std::cout << "Num Cent: " << num_centroids << std::endl; 
    std::cout << "result size: " << ((num_data*num_centroids*4)/1048576) << " Mbytes" << std::endl; 
    data2 = new float[dim*num_data]; 
    centroids2 = new float[dim*num_centroids]; 
    dist2 = new float[num_data*num_centroids]; 
    for (int i = 0; i < dim*num_data; i++) data2[i] = rand()/(float)RAND_MAX; 
    for (int i = 0; i < dim*num_centroids; i++) centroids2[i] = rand()/(float)RAND_MAX; 
    unsigned long long dtime = dtime_usec(0); 
    eucl_dist_cpu(centroids2, data2, dist2, num_centroids, dim, num_data, 0); 
    dtime = dtime_usec(dtime); 
    std::cout << "cpu time: " << dtime/(float)USECPSEC << "s" << std::endl; 
    thrust::host_vector<float> h_data2(data2, data2 + (dim*num_data)); 
    thrust::host_vector<float> h_centr2(centroids2, centroids2 + (dim*num_centroids)); 
    thrust::host_vector<float> h_dist2(num_data*num_centroids); 
    dtime = dtime_usec(0); 
    unsigned long long ctime = eucl_dist_thrust(h_centr2, h_data2, h_dist2, num_centroids, dim, num_data, 0); 
    dtime = dtime_usec(dtime); 
    std::cout << "gpu total time: " << dtime/(float)USECPSEC << "s, gpu compute time: " << ctime/(float)USECPSEC << "s" << std::endl; 
    if (!verify(dist2, &(h_dist2[0]), num_data*num_centroids)) {std::cout << "Verification failure." << std::endl; return 1;} 
    std::cout << "Success!" << std::endl; 

    return 0; 

} 

注：

1. 的代碼被設置爲做2次，短一個使用數據設置類似於你的，與打印進行視覺檢查。然後可以通過命令行大小參數（質心數量，然後是個人數量）輸入更大的數據集，以便進行基準比較和結果驗證。

2. 相反的是我在評論中說，推力代碼只運行比天真的單線程CPU代碼快約25％。你的旅費可能會改變。

3. 這只是一種方式來思考處理它。我有其他的想法，但沒有足夠的時間充實它們。

4. 的數據集可以變得相當大。現在的代碼旨在限制於dimension*number_of_centroids*number_of_individuals的產品少於20億的數據集。但是，即使你接近這個數字，你也需要一個GPU和CPU，它們都有幾GB的內存。我簡要探討了更大的數據集大小。在不同的地方需要改變一些代碼以從例如int到unsigned long long等等。但是我還沒有提供，因爲我仍在調查該代碼的問題。另外，在GPU上計算歐幾里得距離的非推力相關看，你可能會對this question感興趣。如果你遵循那裏所做的優化順序，它可能會揭示這個推力代碼如何被改進，或者如何使用另一個非推力實現。

5. 對不起，我沒能擠出更多的表現出來。