algorithm - 在 for 循环中向量化查找函数

Question

我有以下代码，它输出 array1 的值小于或等于 array2 的每个元素。两个数组的长度不一样。这个 for 循环非常慢，因为数组很大(~500,000 元素)。仅供引用，两个数组始终按升序排列。

任何帮助使它成为向量操作并加快它的速度的人都将不胜感激。

我正在考虑使用“最近”选项的 interp1() 的某种多步骤过程。然后找到相应的 outArray 大于 array2 的位置，然后以某种方式修复点......但我认为必须有更好的方法。

array2 = [5 6 18 25];
array1 = [1 5 9 15 22 24 31];
outArray = nan(size(array2));
for a =1:numel(array2)
    outArray(a) = array1(find(array1 <= array2(a),1,'last'));
end

返回:

outArray =    
     5     5    15    24

Answer 1

这是一种可能的向量化:

[~,idx] = max(cumsum(bsxfun(@le, array1', array2)));
outArray = array1(idx);

---

编辑:

在最近的版本中，得益于 JIT 编译，MATLAB 在执行良好的旧式非矢量化循环方面变得相当出色。

下面是一些类似于您的代码，它利用了两个数组已排序的事实(因此，如果 pos(a) = find(array1<=array2(a), 1, 'last') 那么我们保证在下一次迭代中计算的 pos(a+1) 将不小于之前的 pos(a) )

pos = 1;
idx = zeros(size(array2));
for a=1:numel(array2)
    while pos <= numel(array1) && array1(pos) <= array2(a)
        pos = pos + 1;
    end
    idx(a) = pos-1;
end
%idx(idx==0) = [];      %# in case min(array2) < min(array1)
outArray = array1(idx);

注意:注释行处理 array2 的最小值小于 array1 的最小值(即 find(array1<=array2(a)) 为空时)的情况

我对目前发布的所有方法进行了比较，这确实是最快的方法。长度为 N=5000 的向量的计时(使用 TIMEIT 函数执行)为:

0.097398     # your code
0.39127      # my first vectorized code
0.00043361   # my new code above
0.0016276    # Mohsen Nosratinia's code

这里是 N=500000 的时间:

(? too-long) # your code
(out-of-mem) # my first vectorized code
0.051197     # my new code above
0.25206      # Mohsen Nosratinia's code

.. 从您报告的最初 10 分钟减少到 0.05 秒，这是一个相当不错的改进!

如果你想重现结果，这里是测试代码:

function [t,v] = test_array_find()
    %array2 = [5 6 18 25];
    %array1 = [1 5 9 15 22 24 31];
    N = 5000;
    array1 = sort(randi([100 1e6], [1 N]));
    array2 = sort(randi([min(array1) 1e6], [1 N]));

    f = {...
        @() func1(array1,array2);   %# Aero Engy
        @() func2(array1,array2);   %# Amro
        @() func3(array1,array2);   %# Amro
        @() func4(array1,array2);   %# Mohsen Nosratinia
    };

    t = cellfun(@timeit, f);
    v = cellfun(@feval, f, 'UniformOutput',false);
    assert( isequal(v{:}) )
end

function outArray = func1(array1,array2)
    %idx = arrayfun(@(a) find(array1<=a, 1, 'last'), array2);
    idx = zeros(size(array2));
    for a=1:numel(array2)
        idx(a) = find(array1 <= array2(a), 1, 'last');
    end
    outArray = array1(idx);
end

function outArray = func2(array1,array2)
    [~,idx] = max(cumsum(bsxfun(@le, array1', array2)));
    outArray = array1(idx);
end

function outArray = func3(array1,array2)
    pos = 1;
    lastPos = numel(array1);
    idx = zeros(size(array2));
    for a=1:numel(array2)
        while pos <= lastPos && array1(pos) <= array2(a)
            pos = pos + 1;
        end
        idx(a) = pos-1;
    end
    %idx(idx==0) = [];      %# in case min(array2) < min(array1)
    outArray = array1(idx);
end

function outArray = func4(array1,array2)
    [~,I] = sort([array1 array2]);
    a1size = numel(array1);
    J = find(I>a1size);
    outArray = nan(size(array2));
    for k=1:numel(J),
        if  I(J(k)-1)<=a1size,
            outArray(k) = array1(I(J(k)-1));
        else
            outArray(k) = outArray(k-1);
        end
    end
end