performance - 在 MATLAB 中用 bsxfun 替换 repmat

在以下功能中，我想进行一些更改以使其快速。它本身很快，但我必须在 for 循环中多次使用它，所以它需要很长时间。我想如果我用 bsxfun 替换 repmat 会使其更快，但我不确定。我该如何进行这些替换

function out = lagcal(y1,y1k,source)
kn1 = y1(:);
kt1 = y1k(:);

kt1x = repmat(kt1,1,length(kt1));

eq11 = 1./(prod(kt1x-kt1x'+eye(length(kt1))));
eq1 = eq11'*eq11;

dist = repmat(kn1,1,length(kt1))-repmat(kt1',length(kn1),1);
[fixi,fixj] = find(dist==0); dist(fixi,fixj)=eps;
mult = 1./(dist);

eq2 = prod(dist,2);
eq22 = repmat(eq2,1,length(kt1));
eq222 = eq22 .* mult;

out = eq1 .* (eq222'*source*eq222);
end

它真的加快了我的功能吗？

最佳答案

介绍和代码更改
函数代码中使用的所有 repmat 用法都是将输入扩展到大小，以便稍后可以执行涉及这些输入的数学运算。这是为 bsxfun 量身定制的情况。可悲的是，虽然函数代码的真正瓶颈似乎是别的东西。继续讨论代码的所有性能相关方面。repmat 替换为 bsxfun 的代码如下所示，替换后的代码
保留作为比较的评论 -

function out = lagcal(y1,y1k,source)

kn1 = y1(:);
kt1 = y1k(:);

%//kt1x = repmat(kt1,1,length(kt1));
%//eq11 = 1./(prod(kt1x-kt1x'+eye(length(kt1)))) %//'
eq11 = 1./prod(bsxfun(@minus,kt1,kt1.') + eye(numel(kt1))) %//'

eq1 = eq11'*eq11; %//'

%//dist = repmat(kn1,1,length(kt1))-repmat(kt1',length(kn1),1) %//'
dist = bsxfun(@minus,kn1,kt1.') %//'

[fixi,fixj] = find(dist==0);

dist(fixi,fixj)=eps;
mult = 1./(dist);

eq2 = prod(dist,2);

%//eq22 = repmat(eq2,1,length(kt1));
%//eq222 = eq22 .* mult
eq222 = bsxfun(@times,eq2,mult)

out = eq1 .* (eq222'*source*eq222); %//'

return; %// Better this way to end a function

可以在此处添加更多修改。在最后一行，我们可以做
如下所示，但时序结果并没有显示出巨大的好处
用它 -

out = bsxfun(@times,eq11.',bsxfun(@times,eq11,eq222'*source*eq222))

这将避免在原始代码中较早完成的 eq1 计算，因此您可以通过这种方式节省更多时间。
基准测试
对 bsxfun 修改后的代码部分与原始代码进行基准测试
接下来将讨论基于 repmat 的代码。
基准测试代码

N_arr = [50 100 200 500 1000 2000 3000]; %// array elements for N (datasize)
blocks = 3;
timeall = zeros(2,numel(N_arr),blocks);

for k1 = 1:numel(N_arr)
    N = N_arr(k1);
    y1 = rand(N,1);
    y1k = rand(N,1);
    source = rand(N);

    kn1 = y1(:);
    kt1 = y1k(:);

    %% Block 1 ----------------
    block = 1;
    f = @() block1_org(kt1);
    timeall(1,k1,block) = timeit(f);
    clear f

    f = @() block1_mod(kt1);
    timeall(2,k1,block) = timeit(f);
    eq11 = feval(f);
    clear f
    %% Block 1 ----------------

    eq1 = eq11'*eq11; %//'

    %% Block 2 ----------------
    block = 2;
    f = @() block2_org(kn1,kt1);
    timeall(1,k1,block) = timeit(f);
    clear f

    f = @() block2_mod(kn1,kt1);
    timeall(2,k1,block) = timeit(f);
    dist = feval(f);
    clear f
    %% Block 2 ----------------

    [fixi,fixj] = find(dist==0);

    dist(fixi,fixj)=eps;
    mult = 1./(dist);

    eq2 = prod(dist,2);

    %% Block 3 ----------------
    block = 3;
    f = @() block3_org(eq2,mult,length(kt1));
    timeall(1,k1,block) = timeit(f);
    clear f

    f = @() block3_mod(eq2,mult);
    timeall(2,k1,block) = timeit(f);
    clear f
    %% Block 3 ----------------

end

%// Display benchmark results
figure,
for k2 = 1:blocks
    subplot(blocks,1,k2),
    title(strcat('Block',num2str(k2),' results :'),'fontweight','bold'),hold on
    plot(N_arr,timeall(1,:,k2),'-ro')
    plot(N_arr,timeall(2,:,k2),'-kx')
    legend('REPMAT Method','BSXFUN Method')
    xlabel('Datasize (N) ->'),ylabel('Time(sec) ->')
end

关联函数

function out = block1_org(kt1)
kt1x = repmat(kt1,1,length(kt1));
out = 1./(prod(kt1x-kt1x'+eye(length(kt1))));
return;

function out = block1_mod(kt1)
out = 1./prod(bsxfun(@minus,kt1,kt1.') + eye(numel(kt1)));
return;

function out = block2_org(kn1,kt1)
out = repmat(kn1,1,length(kt1))-repmat(kt1',length(kn1),1);
return;

function out = block2_mod(kn1,kt1)
out = bsxfun(@minus,kn1,kt1.');
return;

function out = block3_org(eq2,mult,length_kt1)
eq22 = repmat(eq2,1,length_kt1);
out = eq22 .* mult;
return;

function out = block3_mod(eq2,mult)
out = bsxfun(@times,eq2,mult);
return;

结果

结论
基于 bsxfun 的代码显示了 2x 比基于 repmat 的代码加速，这是令人鼓舞的。但是对不同数据大小的原始代码的分析显示，最后一行中的多个矩阵乘法似乎占用了函数代码的大部分运行时间，这在 MATLAB 中被认为非常有效。除非您有办法通过使用其他数学技术来避免这些乘法运算，否则它们看起来像是瓶颈。

关于performance - 在 MATLAB 中用 bsxfun 替换 repmat，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/24921424/