通过在 MATLAB 中对 32 位整数进行运算,向量化 hashing/ranking 固定大小的整数组合

Vectorized hashing/ranking of integer combinations of fixed size via operations on 32-bit integers in MATLAB

我在 MATLAB 中动态创建了巨大的 tables/matrices 不同的第一维,其行表示(排序的)6 阶 1-50 范围内的整数组合。

我想为每个组合分配一个唯一值(哈希、排名),以便我可以检查相同的组合是否出现在不同的表中。不同的组合不允许分配相同的值,即没有冲突。我必须在很多这样的表之间做很多这样的比较。因此,出于性能原因,我想通过矢量化 uint32 操作来实现这一点,使其适用于 MATLAB 中的 GPU 加速。

到目前为止我想到的事情:

  1. 词典排序:不知道如何很好地向量化标准的快速递归算法,唯一的选择似乎是 parfor 它通过行,这比其他选项慢。 IIRC,直接显式公式,虽然可向量化,但需要计算二项式,这又需要 log Gamma 函数以避免巨大的阶乘 + double 类型以避免碰撞,如果我没记错的话,即更慢因为它是 'very numerical'.
  2. Cantor 配对函数:可以连续应用 Cantor 配对,这很好,因为它是一个多项式表达式,但它产生的数字远远超过 uint32,并且肯定比其他选项慢。
  3. Base 51(没有双关语意)整数:发送一个 combination/row 向量 (x_1,...,x_6)x_1 + x_2 * 51 + ... + x_6 * 51^5。这是我目前最快的。它很容易向量化,但不幸的是,对于 50 个元素的 rank-6 组合,仍然需要 uint64double,这比 uint32single 类型的操作要慢。

So, I guess, I am looking for a 'clever' injective function on these combinations that computes within the uint32 range and is also well vectorizable (in MATLAB).

如有任何帮助,我们将不胜感激!

编辑: 这是在 uint32singledouble 中对排名和搜索进行基准测试的例程。我使用 MATLAB 的 gputimeit 来生成准确的结果。

% testing parameters:
N = 26;
ord = 5;

% base-(N+1):
basevec_uint32 = gpuArray(repmat(uint32(N+1),1,ord).^cast(0:ord-1,'uint32'));
basevec_single = cast(basevec_uint32,'single').';
basevec_double = cast(basevec_uint32,'double').';

% highest hash value:
max_hash_value = gpuArray(cast(N-ord+1:N,'single'))*basevec_single

% benchmark GPU-accelerated base-(N+1) ranking for uint32, single, and double:
X_uint16 = randi(N,15000000,ord,'uint16','gpuArray');
X_uint32 = cast(X_uint16,'uint32');
X_single = cast(X_uint16,'single');
X_double = cast(X_uint16,'double');

Y_uint16 = randi(N,5000000,ord,'uint16','gpuArray');
Y_uint32 = cast(Y_uint16,'uint32');
Y_single = cast(Y_uint16,'single');
Y_double = cast(Y_uint16,'double');

ranking_uint32 = @() sum(bsxfun(@times,X_uint32,basevec_uint32),2,'native');
ranking_single = @() X_single*basevec_single;
ranking_double = @() X_double*basevec_double;

disp('ranking in uint32:'); gputimeit(ranking_uint32,1)
disp('ranking in single:'); gputimeit(ranking_single,1)
disp('ranking in double:'); gputimeit(ranking_double,1)

% benchmark GPU-accelerated searching in uint32, single, and double matrices:
X_uint32_ranks = myfun_uint32(); Y_uint32_ranks = sum(bsxfun(@times,Y_uint32,basevec_uint32),2,'native');
X_single_ranks = myfun_single(); Y_single_ranks = Y_single*basevec_single;
X_double_hash = myfun_double(); Y_double_ranks = Y_double*basevec_double;

search_uint32 = @() ismember(X_uint32_ranks,Y_uint32_ranks);
search_single = @() ismember(X_single_ranks,Y_single_ranks);
search_double = @() ismember(X_double_ranks,Y_double_ranks);

disp('searching in uint32:'); gputimeit(search_uint32,1)
disp('searching in single:'); gputimeit(search_single,1)
disp('searching in double:'); gputimeit(search_double,1)

我使用的是 nVidia GTX 1060,它在 Windows 下给出了以下基准测试结果(不过,我不知道 WDDM 对 gputimeit 的影响有多大):

>> test1

max_rank_value =

gpuArray single

14327680

ranking in uint32:

ans =

0.0129

ranking in single:

ans =

0.0063

ranking in double:

ans =

0.0108

searching in uint32:

ans =

0.1572

searching in single:

ans =

0.1577

searching in double:

ans =

0.1599

Conclusion: The choice of type does not impact the GPU-accelerated search via ismember (expected), however there is a noticeable difference in the ranking speed coming from the matrix multiplication. Ranking in single is almost twice as fast as ranking in double, while ranking in uint32 is nearly the same as ranking in double. MATLAB does not currently support GPU-accelerated uint32 matrix multiplication, so I am using a substitute function for this. Thus, it is much more reasonable to try to encode combinations of 50 elements of order 6 as single values, which, if I am not mistaken, has just enough bits to accomodate all 15890700 combinations as single integers.

你的最后一个想法几乎已经足够了,所以你只需要根据顺序挤出一些来让它超过标准。因为整个序列是有序的,所以每一对也是有序的。因此,使用 50×50 查找 table 将排序的 (1st,2nd)、(3rd,4th)、(5th,6th) 对映射到 0-1274 之间的数字。

或者如果您不想要 table,可以使用相当简单的显式函数将 j>=i 的对 (i,j) 映射到线性索引。查看上三角或下三角矩阵索引以获取有关这些的详细信息。 (这将是一些类似的东西 n*(n+1)/2 - (n-i)*(n-i-1)/2 + j 根据 base-0 或 base-1 索引,有一些 +/-1 被抛出,在你的情况下 n=50,但我敢肯定我会把它写错。)

无论如何,一旦你得到三个数字 0-1274,base-1275 的想法将适合 uint32