Can you do this calculation with bsxfun instead of repmat?

Question

Henrik 2014-9-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/156327-can-you-do-this-calculation-with-bsxfun-instead-of-repmat

评论： Guillaume 2014-9-27

Hi there

I am trying to optimize some code, and have come quite far with repmat. I think it could become even faster with bsxfun, but I am not sure how to do it. Below are examples of the first implementation with simple for loops and my repmat solution.

L1 and L2 are typically in the range given here, but can potentially be up to around 10 times larger, usually with L2 being the largest.

Thanks in advance!

Henrik

L1=32; %length of f_1, f_2, e_1 and e_2
L2=50; %length of o
%Create some random test data
e_1=rand(L1,1);
f_1=rand(L1,1);
e_2=rand(L1,1);
f_2=rand(L1,1);
o=rand(L2,1);
eta=1e-3; % this just needs to be a small number
tic
Li1=zeros(L1,L1,L2); %preallocate
for n=1:L1;
  for m=1:L1;
      Li1(n,m,:)=(f_1(n)+f_2(m)-1)./(o+1i*eta-e_1(n)-e_2(m));
  end;
end
slowt=toc;
% Same calculation using repmat:
f_1=reshape(f_1,[L1 1 1]);
f_2=reshape(f_2,[1 L1 1 ]);
e_1=reshape(e_1,[L1 1 1]);
e_2=reshape(e_2,[1 L1 1]);
o=reshape(o,[1 1 L2]);
tic
Li2=(repmat(f_1,[1 L1 L2])+repmat(f_2,[L1 1 L2])-1)./...
  (repmat(o,[L1 L1 1])+1i*eta-repmat(e_1,[1 L1 L2])-repmat(e_2,[L1 1 L2]));
fastt=toc;
slowt
fastt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Guillaume 2014-9-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/156327-can-you-do-this-calculation-with-bsxfun-instead-of-repmat#answer_153020

在 MATLAB Online 中打开

I'm not sure bsxfun would help much as because of the division, you need to use it twice. However, I would definitively use meshgrid or ndgrid for part of the calculation:

[f_1, f_2] = ndgrid(f_1, f_2);
[e_1, e_2] = ndgrid(e_1, e_2);
fsum = f_1 + f_2 - 1;
esum = e_1 + e_2

To use bsxfun, you first need to modify o so that the values are in the third dimension

o = shiftdim(o, -2);

Then you can use bsxfun to calculate the denominator and again to perform the division. Each time it expands the 2d matrix (fsum or esum) to 3d along o, and o to 3d along the first two dimensions:

denum = bsxfun(@minus, o+1i*eta, esum);
Li3 = bsxfun(@rdivide, fsum, denum);

3 个评论
显示 1更早的评论隐藏 1更早的评论

Guillaume 2014-9-27

在 MATLAB Online 中打开

Actually, I was so focussed on using bsxfun that I missed the forest for the tree. You don't want to use bsxfun, just ndgrid:

[f_1, f_2] = ndgrid(f_1, f_2, o); %don't want the third output yet
[e_1, e_2, o] = ndgrid(e_1, e_2, o); %want it now
Li3 = (f_1 + f_2 - 1) ./ (o + 1i*eta - e_1 - e_2);

It's pretty much equivalent to your repmat, in fewer lines.

Guillaume 2014-9-27

在 MATLAB Online 中打开

And the way to do it with just bsxfun would be:

f_2 = shiftdim(f_2, -1);
e_2 = shiftdim(e_2, -1);
o = shiftdim(o, -2);
num = bsxfun(@plus, f_1, f_2) + 1;
denum = bsxfun(@minus, o + 1i*eta, bsxfun(@plus, e_1, e_2));
Li2 = bsxfun(@rdivide, num, denum);

A bit convoluted.

请先登录，再进行评论。