matlab中双for循环的矢量化求和
matlab中双for循环的矢量化求和
提问于 2016-09-06 08:39:59
对于非常大的数组(10k X 10k)或更多，我有两个for循环。显然，这个程序的一部分是一个巨大的瓶颈和非常耗时的任务。
共有4个阵列：vm(10000,1)、va(10000,1)、yr(10000,10000)和yi(10000,10000)
for i = 1: 10000
psum = 0;
for j = 1: 10000
psum = psum + vm(i)*vm(j)*(yr(i,j)*cos(va(i)-va(j)) + yi(i,j)*sin(va(i)-va(j)));
end
pcal(i) = psum;
end回答 3
Stack Overflow用户
回答已采纳
发布于 2016-09-06 15:53:54
在您的例子中，一次计算总和是很简单的。基本上，您创建的数组中的元素分别是vm和va的适当乘积和差(使用bsxfun)，然后是各行的逐个元素乘法和求和。
pcal = sum(bsxfun(@times,vm,vm') .* (...
yr.*cos(bsxfun(@minus,va,va')) + ...
yi.*sin(bsxfun(@minus,va,va'))),2);请注意，在矢量化中，您倾向于在内存与CPU周期之间进行权衡。如果您没有足够的RAM，您可能会结束分页，这将使矢量化的解决方案变得缓慢。
Stack Overflow用户
发布于 2016-09-06 19:21:24
最快的选择是@Joans answer。然而，如果你遇到内存问题，这里有一个半矢量化的选项(只有一个循环)：
pcal = zeros(N,1); % N is the 10000 in your example
for m = 1: N
va_va = va(m)-va(1:N);
pcal(m) = sum(vm(m)*vm(1:N).*(yr(m,1:N).'.*cos(va_va)+yi(m,1:N).'.*sin(va_va)));
end下面是该方法的基准测试，以及您的方法和@Joans方法，以及另一个使用ndgrid的方法
function sum_time
N = 10000;
vm = rand(N,1);
va = rand(N,1);
yr = rand(N);
yi = rand(N);
loop_time = timeit(@() loop(N,vm,va,yr,yi))
loop2_time = timeit(@() loop2(N,vm,va,yr,yi))
bsx_time = timeit(@() bsx(vm,va,yr,yi))
ndg_time = timeit(@() ndg(N,vm,va,yr,yi))
end
function pcal = loop(N,vm,va,yr,yi)
pcal = zeros(N,1);
for m = 1: N
psum = 0;
for n = 1: N
psum = psum + vm(m)*vm(n)*(yr(m,n)*cos(va(m)-va(n)) +...
yi(m,n)*sin(va(m)-va(n)));
end
pcal(m) = psum;
end
end
function pcal = loop2(N,vm,va,yr,yi)
pcal = zeros(N,1);
for m = 1: N
va_va = va(m)-va(1:N); % to avoid calculating twice
pcal(m) = sum(vm(m)*vm(1:N).*(yr(m,1:N).'.*cos(va_va)+yi(m,1:N).'.*sin(va_va)));
end
end
function pcal = bsx(vm,va,yr,yi)
pcal = sum(bsxfun(@times,vm,vm') .* (...
yr.*cos(bsxfun(@minus,va,va')) + ...
yi.*sin(bsxfun(@minus,va,va'))),2);
end
function pcal = ndg(N,vm,va,yr,yi)
[n,m] = ndgrid((1:N).',1:N);
yr_t = yr.';
yi_t = yi.';
va_va = va(m(:))-va(n(:));
vmt = vm(m(:)).*vm(n(:));
psum = vmt.*(yr_t(1:N^2).'.*cos(va_va)+yi_t(1:N^2).'.*sin(va_va));
pcal = sum(reshape(psum,N,N)).';
end以及结果(针对N = 10000)：
loop_time =
7.0296
loop2_time =
3.3722
bsx_time =
1.2716
ndg_time =
6.3568因此，减少一次循环可节省约50%的时间。
Stack Overflow用户
发布于 2016-09-07 06:01:36
你可以根据这些trigonometric identities重新表述你的方程。
sin(a-b) = sin a cos b - cos a sin b;
cos(a-b) = cos a cos b + sin a sin b;因此，预先计算正弦和余弦，并在loop或bsxfun.中使用它们。这是循环版本：
yr = rand(10000);
yi = rand(10000);
va = rand(1,10000);
vm = rand(1,10000);
sin_va = sin(va);
cos_va = cos(va);
for i = 1: 10000
pcal(i) = sum(vm(i)*vm.*(yr(i,:).*(cos_va(i) * cos_va + sin_va(i) * sin_va) + yi(i,:).*(sin_va(i) * cos_va - cos_va(i) * sin_va)));
end页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/39339149
复制相关文章
相似问题
腾讯云开发者
