matlab - 如何矢量化在较大矩阵的子集上运行函数的代码？

Question

假设我有以下 9 x 5 矩阵:

myArray = [
   54.7    8.1   81.7   55.0   22.5
   29.6   92.9   79.4   62.2   17.0
   74.4   77.5   64.4   58.7   22.7
   18.8   48.6   37.8   20.7   43.5
   68.6   43.5   81.1   30.1   31.1
   18.3   44.6   53.2   47.0   92.3
   36.8   30.6   35.0   23.0   43.0
   62.5   50.8   93.9   84.4   18.4
   78.0   51.0   87.5   19.4   90.4
];

我有这个矩阵的 11 个“子集”，我需要对每个子集运行一个函数(比如说 max)。可以使用以下逻辑矩阵来标识子集(按列标识，而不是按行标识):

myLogicals = logical([
    0 1 0 1 1
    1 1 0 1 1
    1 1 0 0 0
    0 1 0 1 1
    1 0 1 1 1
    1 1 1 1 0
    0 1 1 0 1
    1 1 0 0 1
    1 1 0 0 1
]);

或通过线性索引:

starts = [2 5 8 10 15 23 28 31 37 40 43]; #%index start of each subset
ends =   [3 6 9 13 18 25 29 33 38 41 45]; #%index end of each subset

第一个子集是 2:3，第二个子集是 5:6，依此类推。

我可以找到每个子集的 max 并将其存储在一个向量中，如下所示:

finalAnswers = NaN(11,1); 
for n=1:length(starts) #%i.e. 1 through the number of subsets
    finalAnswers(n) = max(myArray(starts(n):ends(n)));
end

循环运行后，finalAnswers 包含每个数据子集的最大值:

74.4  68.6  78.0  92.9  51.0  81.1  62.2  47.0  22.5  43.5  90.4

是否可以在不使用 for 循环的情况下获得相同的结果？换句话说，这段代码可以向量化吗？这种方法会比目前的方法更有效吗？

---

编辑:我对建议的解决方案进行了一些测试。我使用的数据是一个 1,510 x 2,185 矩阵，包含 10,103 个子集，长度从 2 到 916 不等，子集长度的标准差为 101.92。

我将每个解决方案包装在 tic;for k=1:1000 [code here] end; toc; 结果如下:

• for 循环方法 --- Elapsed time is 16.237400 seconds.
• Shai 的方法 --- 耗时是 153.707076 秒。
• Dan 的方法 --- 耗时是 44.774121 秒。
• Divakar 的方法 #2 --- 耗时是 127.621515 秒。

注意事项:

• 我还尝试通过将 k=1:1000 for 循环包装在 accumarray 行周围来对 Dan 的方法进行基准测试(因为其余部分可能是理论上只运行一次)。在这种情况下，时间是 28.29秒。
• 对 Shai 的方法进行基准测试，同时保留 lb = ... 行k 循环的时间是 113.48 秒。
• 当我运行 Divakar 的代码时，我得到了 两个的非单一维度
输入数组必须彼此匹配。 bsxfun 行的错误。我通过使用共轭转置(撇号trade_starts(1:starts_extent) 上的运算符 ') 和intv(1:starts_extent) 在调用 bsxfun 的代码行中。我是不确定为什么会出现此错误...

我不确定我的基准测试设置是否正确，但看起来 for 循环在这种情况下实际上运行最快。

Answer 1

一种方法是使用accumarray。不幸的是，为了做到这一点，我们首先需要“标记”您的逻辑矩阵。如果您没有图像处理工具箱，这是一种复杂的方法:

sz=size(myLogicals);
s_ind(sz(1),sz(2))=0;
%// OR: s_ind = zeros(size(myLogicals))

s_ind(starts) = 1;
labelled = cumsum(s_ind(:)).*myLogicals(:);        

所以这就是 Shai 的 bwlabeln 实现所做的(但这将是 1-by-numel(myLogicals) 的形状，而不是到 size(myLogicals) 的形状)

现在你可以使用accumarray了:

accumarray(labelled(myLogicals), myArray(myLogicals), [], @max)

否则尝试可能会更快

result = accumarray(labelled+1, myArray(:), [], @max);
result = result(2:end)

这是完全矢量化的，但值得吗？您必须针对循环解决方案进行速度测试才能知道。