在matlab中用双for循环矢量化求和答案

【问题标题】：vectorizing summation with double for loops in matlab在matlab中用双for循环矢量化求和
【发布时间】：2016-09-07 06:01:36
【问题描述】：

对于非常大的数组（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

【问题讨论】：

尝试从矢量化here 和逐元素乘法here 中获得一些灵感，并尝试为您的问题提供自己的第一个解决方案，您将从那里获得帮助。您可以使用 tic-toc 模式在 matlab 中测量时间
我还建议不要直接使用值 10000，而是使用 size 或 numel。如果明天你的输入向量有 5000 个元素，你会得到一个错误，但如果向量有 20000 个元素，你几乎不会注意到其中一半没有被计算。

标签： matlab vectorization

【解决方案1】：

在您的情况下，一次性计算总和很简单。基本上，您创建数组，其中元素分别是 vm 和 va 的适当乘积和差异（使用 bsxfun），然后是逐行乘法和求和。

pcal = sum(bsxfun(@times,vm,vm') .* (...
    yr.*cos(bsxfun(@minus,va,va')) + ...
    yi.*sin(bsxfun(@minus,va,va'))),2);

请注意，在 vecorization 中，您倾向于交换内存与 CPU 周期。如果您没有足够的 RAM，则可能会出现分页，这会减慢矢量化解决方案的爬行速度。

【讨论】：

非常感谢 Jonas，这段代码运行良好。这是我的问题的答案

【解决方案2】：

最快的选择是@Joans 回答。但是，如果您遇到内存问题，这里有一个半向量化选项（只有一个循环）：

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% 的时间。

【讨论】：

感谢 EBH，感谢您提供详细选项的解释

【解决方案3】：

您可以根据论文trigonometric identities 重新制定您的方程式：

sin(a-b) = sin a cos b - cos a sin b;
cos(a-b) = cos a cos b + sin a sin b;

所以预先计算正弦和余弦并在循环或 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

【讨论】：