矩阵的矢量化列逻辑索引

Question

我正在MATLAB中对矩阵进行逐列逻辑索引。例如：

tic
N = 5*10^6;
input = randi(100,N,12);
output = zeros(N,1);
d_sn2 = randi(25,N,12);
d_sd2 = randi(25,N,12);
LL1 = randi(8,N,12);
UL1 = randi([12,20],N,12);
LL2 = randi(8,N,12);
UL2 = randi([12,20],N,12);
for p = 1:N
temp = zeros(12,12);
for i = 1:12
I2 = (d_sn2(:,i)>LL1(p,i) & d_sn2(:,i)<UL1(p,i)) & (d_sd2(:,i)>LL2(p,i) & d_sd2(:,i)<UL2(p,i));
temp(i,:) = mean(input(I2,:));
end
output(p) = max(temp(:));
end
toc               

我想知道我是否可以将此操作矢量化或使其更快？

Answer 1

内循环中I2的计算可以很容易地矢量化。此：

temp = zeros(12,12);
for i = 1:12
I2 = (d_sn2(:,i)>LL1(p,i) & d_sn2(:,i)<UL1(p,i)) & (d_sd2(:,i)>LL2(p,i) & d_sd2(:,i)<UL2(p,i));
temp(i,:) = mean(input(I2,:));
end

与此相同：

I2 = d_sn2>LL1(p,:) & d_sn2<UL1(p,:) & d_sd2>LL2(p,:) & d_sd2<UL2(p,:);
temp = zeros(12,12);
for i = 1:12
temp(i,:) = mean(input(I2(:,i),:));
end

此代码使用隐式单例扩展，如果您有R2016b之前的MATLAB版本，则需要使用bsxfun:bsxfun(@gt,d_sn2,LL1(p,:))等编写每个>(gt(和<(lt(调用。

不幸的是，对input进行索引要向量化困难得多。因为i的每次迭代都访问不同数量的input元素，所以没有简单的方法可以创建没有循环的矩阵temp。我尝试过的几种方法都比循环代码慢得多。

如果你使用的是一个相当新的MATLAB版本，你的代码将非常有效。MATLAB的解释器使用JIT(实时编译器(，使循环速度不再像以前那么慢。例如，一个平凡的循环添加矩阵所有元素的速度只比使用函数sum慢2-3倍。在过去，这可能要慢100倍。因此，矢量化的好处与过去不一样了。再加上你正在使用的非常大的数组，这意味着矢量化将是一种令人讨厌的方法，因为矢量化通常意味着创建更大的中间矩阵。