迭代模拟存在严重的性能问题

Question

最近在实现一个模拟算法时，我遇到了性能问题。我成功地找到了瓶颈函数（显然，它是内部调用 arrayfun 的速度变慢）：

function sim = simulate_frequency(the_f,k,n)

    r = rand(1,n); % 
    x = arrayfun(@(x) find(x <= the_f,1,'first'),r);
    sim = (histcounts(x,[1:k Inf]) ./ n).';

end

它在代码的其他部分中被用作如下:

h0 = zeros(1,sims);

for i = 1:sims
    p = simulate_frequency(the_f,k,n);
    h0(i) = max(abs(p - the_p));
end

以下是一些可能的取值：

% Test Case 1
sims = 10000;
the_f = [0.3010; 0.4771; 0.6021; 0.6990; 0.7782; 0.8451; 0.9031; 0.9542; 1.0000];
k = 9;
n = 95;

% Test Case 2
sims = 10000;
the_f = [0.0413; 0.0791; 0.1139; 0.1461; 0.1760; 0.2041; 0.2304; 0.2552; 0.2787; 0.3010; 0.3222; 0.3424; 0.3617; 0.3802; 0.3979; 0.4149; 0.4313; 0.4471; 0.4623; 0.4771; 0.4913; 0.5051; 0.5185; 0.5314; 0.5440; 0.5563; 0.5682; 0.5797; 0.5910; 0.6020; 0.6127; 0.6232; 0.6334; 0.6434; 0.6532; 0.6627; 0.6720; 0.6812; 0.6901; 0.6989; 0.7075; 0.7160; 0.7242; 0.7323; 0.7403; 0.7481; 0.7558; 0.7634; 0.7708; 0.7781; 0.7853; 0.7923; 0.7993; 0.8061; 0.8129; 0.8195; 0.8260; 0.8325; 0.8388; 0.8450; 0.8512; 0.8573; 0.8633; 0.8692; 0.8750; 0.8808; 0.8864; 0.8920; 0.8976; 0.9030; 0.9084; 0.9138; 0.9190; 0.9242; 0.9294; 0.9344; 0.9395; 0.9444; 0.9493; 0.9542; 0.9590; 0.9637; 0.9684; 0.9731; 0.9777; 0.9822; 0.9867; 0.9912; 0.9956; 1.000];
k = 90;
n = 95;

标量sims必须在1000和1000000之间。累积频率向量the_f中不超过100个元素。标量k表示the_f中的元素数量。最后，标量n表示经验样本向量中的元素数量，甚至可以非常大（据我所知，高达10000个元素）。

有什么方法可以改善此过程的计算时间吗？

- Tommaso Belluzzo

Answer 1

在第二个测试用例中，这对我来说似乎稍微快了一点，而不是第一个。对于更长的the_f和更大的n值，时间差异可能会更大。

function sim = simulate_frequency(the_f,k,n)
    r = rand(1,n); % 
    [row,col] = find(r <= the_f); % Implicit singleton expansion going on here!
    [~,ind] = unique(col,'first');
    x = row(ind);
    sim = (histcounts(x,[1:k Inf]) ./ n).';
end

我正在使用r <= the_f的隐式单例展开，如果您使用旧版本的MATLAB（但您知道该怎么做），请使用bsxfun。

find然后返回所有位置的行和列，其中r大于the_f。unique查找结果的索引，以获取每列的第一个元素。

^{来源：MATLAB Answers上的Andrei Bobrov}

---

另一种选项（源自这个答案）更短，但在我看来也更加晦涩：

mask = r <= the_f;
[x,~] = find(mask & (cumsum(mask,1)==1));

- Cris Luengo

---

我喜欢即使你不使用bsxfun :-P，你也添加了“隐式单例扩展正在进行中”的注释。 - Luis Mendo

Answer 2

如果我需要更好的性能，我会避免使用arrayfun。即使是这个for循环也更快：

function sim = simulate_frequency(the_f,k,n)

    r = rand(1,n); % 
    for i = 1:numel(r)
        x(i) = find(r(i)<the_f,1,'first');
    end
    sim = (histcounts(x,[1:k Inf]) ./ n).';

end

使用第一组样本数据运行10000个模拟，得到以下计时结果。

您的arrayfun函数：

>Elapsed time is 2.848206 seconds.

for循环函数：

>Elapsed time is 0.938479 seconds.

受Cris Luengo答案的启发，我建议如下：

function sim = simulate_frequency(the_f,k,n)

    r = rand(1,n); % 
    x = sum(r > the_f)+1;
    sim = (histcounts(x,[1:k Inf]) ./ n)';

end

时间：

>Elapsed time is 0.264146 seconds.

- Anthony

---

sum(r > the_f)+1 -- 真聪明！ - Cris Luengo

Answer 3

您可以使用histcounts函数，并将其输入设置为r:

r = rand(1,n);
sim = (histcounts(r,[-inf ;the_f]) ./ n).';

如果使用histc而不是histcounts，整个模拟可以向量化：

r = rand(n,sims);
p = histc(r, [-inf; the_f],1);
p = [p(1:end-2,:) ;sum(p(end-1:end,:))]./n;
h0 = max(abs(p-the_p(:)));    %h0 = max(abs(bsxfun(@minus,p,the_p(:))));

- rahnema1