用MATLAB中的bsxfun替换repmat [英] replace repmat with bsxfun in MATLAB
问题描述
$ p $
function out = lagcal（y1，y1k，source）
kn1 = y1 ）;
kt1 = y1k（:);
kt1x = repmat（kt1,1，length（kt1））; （b）（b）（b）（b）（b）
eq1 = eq11'* eq11;
$ b $ 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 $ b $ 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）
$ b $ out = eq1。*（eq222 * * source * eq222）; ％//'
return; ％//更好地结束函数
可以在这里添加一个修改。在最后一行，我们可以做
的东西，如下所示，但时间结果并没有显示出巨大的好处
与它 -
out = bsxfun（@ times，eq11。'，bsxfun（@ times，eq11，eq222'* source * eq222））
这样可以避免在原始代码的前面完成 eq1 的计算，所以你可以节省更多的时间 b
标杆管理
bsxfun代码修改部分与原始
repmat接下来讨论基于代码的代码。
标杆代码
N_arr = [50 100 200 500 1000 2000 3000]; ％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
％//显示基准测试结果
figure，
for k2 = 1：blocks
subplot（blocks，1，k2），
title（strcat（'Block'，num2str（k2），' （N_arr，timeall（1，：k2），' - ro'）
plot（N_arr，timeall（2， （'Datasize（N） - >'），ylabel（'Time（sec）'，'k2'，' - kx'）
legend（'REPMAT Method'，'BSXFUN Method'）
xlabel - >'）
end
相关函数 （kt1,1，length（kt1））;这个函数是一个函数，它是一个函数。 （b）（b）（b）（b）
return;
function out = block1_mod（kt1）
out = 1./prod(bsxfun(@minus,kt1,kt1。'）+ eye（numel（kt1）））;
return;
$ b函数out = block2_org（kn1，kt1）
out = repmat（kn1,1，length（kt1）） - repmat（kt1'，length（kn1），1）
return;
函数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;
out = block3_mod（eq2，mult）
out = bsxfun（@ times，eq2，mult）;
return;
结果
结论
bsxfun 基于代码显示 2x 基于repmat的加速这是令人鼓舞的。但是，在不同的数据大小上对原始代码进行剖析，显示在最后一行中的多个矩阵乘法似乎占用了功能代码的大部分运行时间，这在MATLAB中是非常有效的。除非你有办法通过其他数学方法来避免这些乘法运算，否则它们看起来就像是瓶颈。
In the following function i want to make some changes to make it fast. By itself it is fast but i have to use it many times in a for loop so it takes long. I think if i replace the repmat with bsxfun will make it faster but i am not sure. How can i do these replacements
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
Does it really speed up my function?
Introduction and code changes
All the repmat usages used in the function code are to expand inputs to sizes so that later on the mathemtical operations involving these inputs could be performed. This is tailor-made situation for bsxfun. Sadly though the real bottleneck of the function code seems to be something else. Stay on as we discuss all the performance related aspects of the code.
Code with repmat replaced by bsxfun is presented next and the replaced codes
are kept as comments for comparison -
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
One more modification could be added here. In the last line, we could do something like as shown below, but the timing results don't show a huge benefit with it -
out = bsxfun(@times,eq11.',bsxfun(@times,eq11,eq222'*source*eq222))
This would avoid the calculation of eq1 done earlier in the original code, so you would save little more time that way.
Benchmarking
Benchmarking on the bsxfun modified portions of the code versus the original repmat based codes is discussed next.
Benchmarking Code
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
Associated functions
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;
Results
Conclusions
bsxfun based codes show around 2x speedups over repmat based ones which is encouraging. But a profiling of the original code across a varying datasize show the multiple matrix multiplications in the final line seem to be occupying most of the runtime for the function code, which are supposedly very efficient within MATLAB. Unless you have some way to avoid those multiplications by using some other mathematical technique, they look like the bottleneck.
这篇关于用MATLAB中的bsxfun替换repmat的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
