tl;dr: 最快的方法取决于数组的大小。对于小于2
14的数组,下面的方法3 (
accumarray
) 更快。对于大于该大小的数组,下面的方法2 (
histcounts
) 更好。
更新:我还使用了2016b引入的
implicit broadcasting进行了测试,结果几乎与
bsxfun
方法相同,在该方法中没有显着差异(相对于其他方法)。
让我们看看有哪些可用的方法来执行此任务。对于以下示例,我们将假设
X
具有
n
个元素,从1到
n
,我们感兴趣的数组是
M
,它是一个大小可能变化的列数组。我们的结果向量将是
spp
1,其中
spp(k)
是
M
中
k
的数量。尽管我在这里写关于
X
,但在下面的代码中没有明确的实现,我只定义
n = 500
,并且
X
隐式地为
1:500
。
朴素的for
循环
解决此任务最简单和直接的方法是使用一个
for
循环迭代
X
中的元素,并计算等于它的
M
中的元素数量:
function spp = loop(M,n)
spp = zeros(n,1);
for k = 1:size(spp,1);
spp(k) = sum(M==k);
end
end
当然,这样做并不聪明,特别是当只有来自X
的少量元素填充M
时,因此我们最好先寻找那些已经在M
中的元素:
function spp = uloop(M,n)
u = unique(M);
spp = zeros(n,1);
for k = u(u>0).';
spp(k) = sum(M==k);
end
end
通常,在MATLAB中,尽可能利用内置函数是可取的,因为大多数情况下它们更快。我想到了5个选项来实现这一点:
1. 函数
tabulate
函数
tabulate
返回一个非常方便的频率表,乍一看似乎是这个任务的完美解决方案:
function tab = tabi(M)
tab = tabulate(M);
if tab(1)==0
tab(1,:) = [];
end
end
唯一需要修复的是,如果表格中计数了
0
元素(可能
M
中没有零),则需要删除表格的第一行。
2.函数
histcounts
。另一个可以很容易地根据我们的需求进行调整的选项是
histcounts
:
function spp = histci(M,n)
spp = histcounts(M,1:n+1);
end
在这里,为了分别计算1到
n
之间的所有不同元素,我们将边界定义为
1:n+1
,因此
X
中的每个元素都有自己的箱子。我们也可以编写
histcounts(M(M>0),'BinMethod','integers')
,但我已经测试过它,它需要更多时间(尽管它使函数独立于
n
)。
3. 函数
accumarray
接下来我要介绍的选项是使用函数
accumarray
:
function spp = accumi(M)
spp = accumarray(M(M>0),1);
end
在这里,我们将函数
M(M>0)
作为输入,以跳过零,并使用
1
作为
vals
输入来计算所有唯一元素的数量。
我们甚至可以使用二进制运算符
@eq
(即
==
)来查找每种类型的所有元素:
function spp = bsxi(M,n)
spp = bsxfun(@eq,M,1:n);
spp = sum(spp,1);
end
如果我们将第一个输入
M
和第二个
1:n
放在不同的维度中,那么一个是列向量,另一个是行向量,然后函数将比较
M
中的每个元素与
1:n
中的每个元素,并创建一个
length(M)
-by-
n
逻辑矩阵,我们可以对其进行求和以获得所需的结果。
另一种选项类似于
bsxfun
,即使用
ndgrid
函数显式创建所有可能性的两个矩阵:
function spp = gridi(M,n)
[Mx,nx] = ndgrid(M,1:n);
spp = sum(Mx==nx);
end
然后我们将它们进行比较并在列上求和,以得到最终结果。
基准测试
我进行了一项小测试,以找出所有上述方法中最快的方法,我对所有测试定义了 n = 500
。对于某些方法(特别是天真的 for
),n
对执行时间有很大的影响,但这不是问题,因为我们想要针对给定的 n
进行测试。
以下是测试结果:![Timing hist](https://istack.dev59.com/AS61T.webp)
我们可以注意到几件事:
- 有趣的是,最快的方法发生了转变。对于小于214的数组来说,
accumarray
是最快的。对于大于214的数组,histcounts
是最快的。
- 如预期的那样,在两个版本中,朴素的
for
循环都是最慢的,但对于小于28的数组,“unique & for”选项更慢。在大于211 的数组中,ndgrid
成为最慢的,这可能是因为需要在内存中存储非常大的矩阵。
tabulate
在小于29的大小的数组上工作的方法存在一些不规则性。在我进行的所有试验中,这个结果是一致的(模式略有变化)。
(由于在更高的值上会导致计算机卡顿,并且趋势已经非常明显),所以bsxfun
和ndgrid
曲线被截断了)
另外,请注意y轴为log10,因此单位减少(比如对于大小为219的数组,在accumarray
和histcounts
之间),意味着操作速度快了10倍。
欢迎在评论中提出改进此测试的建议,如果您有另一种概念上不同的方法,请建议为答案。
代码
这里是所有函数包装在一个计时函数中:
function out = timing_hist(N,n)
M = randi([0 n],N,1);
func_times = {'for','unique & for','tabulate','histcounts','accumarray','bsxfun','ndgrid';
timeit(@() loop(M,n)),...
timeit(@() uloop(M,n)),...
timeit(@() tabi(M)),...
timeit(@() histci(M,n)),...
timeit(@() accumi(M)),...
timeit(@() bsxi(M,n)),...
timeit(@() gridi(M,n))};
out = cell2mat(func_times(2,:));
end
function spp = loop(M,n)
spp = zeros(n,1);
for k = 1:size(spp,1);
spp(k) = sum(M==k);
end
end
function spp = uloop(M,n)
u = unique(M);
spp = zeros(n,1);
for k = u(u>0).';
spp(k) = sum(M==k);
end
end
function tab = tabi(M)
tab = tabulate(M);
if tab(1)==0
tab(1,:) = [];
end
end
function spp = histci(M,n)
spp = histcounts(M,1:n+1);
end
function spp = accumi(M)
spp = accumarray(M(M>0),1);
end
function spp = bsxi(M,n)
spp = bsxfun(@eq,M,1:n);
spp = sum(spp,1);
end
function spp = gridi(M,n)
[Mx,nx] = ndgrid(M,1:n);
spp = sum(Mx==nx);
end
这是运行此代码并生成图形的脚本:
N = 25;
func_times = zeros(N,5);
for n = 1:N
func_times(n,:) = timing_hist(2^n,500);
end
hold on
mark = 'xo*^dsp';
for k = 1:size(func_times,2)
plot(1:size(func_times,1),log10(func_times(:,k).*1000),['-' mark(k)],...
'MarkerEdgeColor','k','LineWidth',1.5);
end
hold off
xlabel('Log_2(Array size)','FontSize',16)
ylabel('Log_{10}(Execution time) (ms)','FontSize',16)
legend({'for','unique & for','tabulate','histcounts','accumarray','bsxfun','ndgrid'},...
'Location','NorthWest','FontSize',14)
grid on
1 这个奇怪的名称来源于我的领域,生态学。我的模型是元胞自动机,通常模拟虚拟空间中的个体生物(上面的M
)。这些个体属于不同的物种(因此称为spp
),共同组成所谓的“生态群落”。群落的“状态”由每个物种的个体数量组成,即本答案中的spp
向量。在这些模型中,我们首先为个体定义一个物种库(上面的X
),而群落状态考虑了物种库中的所有物种,而不仅仅是存在于M
中的物种。
9
是一个有效的数字,恰巧没有出现呢? - Dev-iL