简介和代码更改
函数代码中使用的所有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
基于 repmat 的代码显示2x
了比基于 repmat 的代码更快的速度,这是令人鼓舞的。但是对不同数据大小的原始代码的分析表明,最后一行中的多个矩阵乘法似乎占据了函数代码的大部分运行时间,这在 MATLAB 中被认为是非常有效的。除非你有办法通过使用其他数学技术来避免这些乘法,否则它们看起来就像是瓶颈。