加速嵌套的for循环

Question

我一直在努力加速以下函数，但一直没有结果：

function beta = beta_c(k,c,gamma)
beta = zeros(size(k));
E = @(x) (1.453*x.^4)./((1 + x.^2).^(17/6));
for ii = 1:size(k,1)
    for jj = 1:size(k,2)
        E_int = integral(E,k(ii,jj),10000);
        beta(ii,jj) = c*gamma/(k(ii,jj)*sqrt(E_int));
    end
end
end

到目前为止，我是这样解决的:

function beta = beta_calc(k,c,gamma)
k_1d = reshape(k,[1,numel(k)]);
E_1d =@(k) 1.453.*k.^4./((1 + k.^2).^(17/6));
E_int = zeros(1,numel(k_1d));
parfor ii = 1:numel(k_1d)
E_int(ii) = quad(E_1d,k_1d(ii),10000);
end
beta_1d = c*gamma./(k_1d.*sqrt(E_int));
beta = reshape(beta_1d,[size(k,1),size(k,2)]);
end

在我看来，它并没有真正提高性能。你对此有什么看法？

你介意给个提示吗？

提前感谢你。

编辑

我将介绍一些涉及我的问题的理论背景。 通常，beta的计算如下所示

因此，在一维k数组的简化情况下，可以计算出E_int

E = 1.453.*k.^4./((1 + k.^2).^(17/6));
E_int = 1.5 - cumtrapz(k,E);

或者，作为替代方案

E_int(1) = 1.5;
for jj = 2:numel(k)
E =@(k) 1.453.*k.^4./((1 + k.^2).^(17/6));
E_int(jj) = E_int(jj - 1) - integral(E,k(jj-1),k(jj));
end

然而，k 目前是一个矩阵 k(size1,size2)。

- fpe

---

你最好的选择是将该函数编译为mex文件 - 运行速度会快很多。 - Floris

@fpe 题外话：在当前版本的MATLAB中，我从未见过arrayfun比for循环更快。 - Eitan T

我不熟悉积分，但它不能接受两个矩阵吗？ - Pavan Yalamanchili

@EitanT：顺便说一下，我想优化上述函数的性能。矢量化可能是答案吗？ - fpe

@Pavan：我猜积分需要标量作为积分限制。 - fpe

Answer 1

我喜欢这个问题。

问题：函数integral只接受标量作为积分限制。因此，矢量化计算E_int的过程变得困难。

提示：从k(ii,jj)到无穷大积分相同的函数似乎存在很多冗余......

解决方案：将k的值从小到大排序，然后用E_sort_int(si) = integral( E, sortedK(si), sortedK(si+1) );来积分，其中sortedK(numel(k)+1) = 10000;。然后完整的E_int = cumsum(E_sort_int);（您只需要“撤销”排序并将其重新调整为k的大小即可）。

- Shai

---

实际上，E_int 作为积分限制应该假定为 |k| 和 Inf；顺便说一句，正如你所争论的那样，我可以用另一种方式计算 E_int，我会将其作为对主要问题的编辑发布（请注意，我将在 k 是一维数组的情况下呈现它）。 - fpe

我开始思考，也许编写一个 C 函数来处理我现在在 Matlab 中调用的这个函数会更容易，至少在涉及到integral时是这样。 - fpe

Answer 2

这里还有另一种方法，即并行化，因为使用 spmd 或 parfor 很容易。不要使用 integral，考虑使用 quad，请参阅此链接获取示例...

- bla

---

我刚刚编辑了我的问题，在beta_calc的表达式中加入了parfor。对于一个案例研究magic(100)，执行时间减少了一半。这是相当不错的结果，但我不知道是否足够好。 - fpe

你的提示似乎是目前最合理的解决方案。谢谢。 - fpe