从两个三维多集数组中更快地找到任意两个对应多集的交集的大小

Question

我有两个uint16三维(GPU)阵列A和B在MATLAB中，其中有相同的第二和第三维空间。例如，size(A,1) = 300 000，size(B,1) = 2000，size(A,2) = size(B,2) = 20和size(A,3) = size(B,3) = 100，给出一个数量级的概念。实际上，size(A,3) = size(B,3)是非常大的，比如说~ 1 000 000，但是这些阵列是在外部存储在沿三维切割的小块中。重点是，沿着三维有一个很长的循环(cfg )。因此，其中的代码需要进一步优化(如果可能的话)。此外，可以假定A和B的值远低于65535，但仍然有数百个不同的值。

对于每个i、j和d，行A(i,:,d)和B(j,:,d) 表示相同大小的多个集合，我需要找到最大的公共子集的大小(多子集？)其中，即它们的交集的大小为多集。此外，可以假定B的行是排序的。

例如，如果[2 3 2 1 4 5 5 5 6 7]和[1 2 2 3 5 5 7 8 9 11]分别是两个这样的多集，那么它们的多集交集是[1 2 2 3 5 5 7]，它的大小为7(7个元素作为一个多集)。

我目前正在使用以下例程来执行此操作：

s = 300000; % 1st dim. of A
n = 2000; % 1st dim. of B
c = 10; % 2nd dim. of A and B
depth = 10; % 3rd dim. of A and B (corresponds to a batch of size 10 of A and B along the 3rd dim.)
N = 100; % upper bound on the possible values of A and B

A = randi(N,s,c,depth,'uint16','gpuArray');
B = randi(N,n,c,depth,'uint16','gpuArray');

Sizes_of_multiset_intersections = zeros(s,n,depth,'uint8'); % too big to fit in GPU memory together with A and B
for d=1:depth
    A_slice = A(:,:,d);
    B_slice = B(:,:,d);
    unique_B_values = permute(unique(B_slice),[3 2 1]); % B is smaller than A

    % compute counts of the unique B-values for each multiset:
    A_values_counts = permute(sum(uint8(A_slice==unique_B_values),2,'native'),[1 3 2]);
    B_values_counts = permute(sum(uint8(B_slice==unique_B_values),2,'native'),[1 3 2]);

    % compute the count of each unique B-value in the intersection:
    Sizes_of_multiset_intersections_tmp = gpuArray.zeros(s,n,'uint8');
    for i=1:n
        Sizes_of_multiset_intersections_tmp(:,i) = sum(min(A_values_counts,B_values_counts(i,:)),2,'native');
    end

    Sizes_of_multiset_intersections(:,:,d) = gather(Sizes_of_multiset_intersections_tmp);
end

我们也可以很容易地修改上面的代码，按照维度3而不是d=1:depth (=批大小1)成批计算结果，但代价是牺牲更大的unique_B_values向量。

由于depth维度很大(即使是在批处理时也是如此)，所以我对外部循环中代码的更快替代方案感兴趣。因此，我的问题是：是否有一种更快(例如更好的矢量化)方法来计算等尺寸多集的交叉口的大小？

Answer 1

免责声明：这不是一个基于GPU的解决方案(不要有一个好的GPU)。我觉得这个结果很有趣，想和大家分享，但是如果你认为应该的话，我可以删除这个答案。

下面是代码的向量化版本，这样就可以摆脱内部循环，而代价是必须处理更大的数组，而这个数组可能太大，无法在内存中使用。

这样做的目的是让矩阵A_values_counts和B_values_counts是三维矩阵，这样调用min(A_values_counts,B_values_counts)就可以一次计算出由于隐式展开而产生的所有数据。在后台，它将创建一个大数组的大小s x n x length(unique_B_values) (可能大多数情况下太大)

为了绕过对大小的约束，在n维数(即B的第一维)上分批计算结果：

tic

nBatches_B = 2000;
sBatches_B = n/nBatches_B;
Sizes_of_multiset_intersections_new = zeros(s,n,depth,'uint8');

for d=1:depth
    A_slice = A(:,:,d);
    B_slice = B(:,:,d);

    % compute counts of the unique B-values for each multiset:    
    unique_B_values = reshape(unique(B_slice),1,1,[]);

    A_values_counts = sum(uint8(A_slice==unique_B_values),2,'native'); % s x 1 x length(uniqueB) array
    B_values_counts = reshape(sum(uint8(B_slice==unique_B_values),2,'native'),1,n,[]); % 1 x n x length(uniqueB) array

    % Not possible to do it all in one go, must split in batches along B

    for ii = 1:nBatches_B

        Sizes_of_multiset_intersections_new(:,((ii-1)*sBatches_B+1):ii*sBatches_B,d) = sum(min(A_values_counts,B_values_counts(:,((ii-1)*sBatches_B+1):ii*sBatches_B,:)),3,'native'); % Vectorized
    
    end

end

toc

下面是一个具有不同批数值的小基准。您可以看到，在400 (批次大小50)的最小值附近，处理时间减少了大约10% (每个点平均超过3次)。(编辑:X轴是批次的数量，而不是批次大小)

​

​

我很想知道它对于GPU阵列的行为是怎样的！