Can you do this calculation with bsxfun instead of repmat?

1 次查看（过去 30 天）
Henrik
Henrik 2014-9-26
评论： 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

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

采纳的回答

Guillaume
Guillaume 2014-9-26
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更早的评论
Guillaume
Guillaume 2014-9-27
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
Guillaume 2014-9-27
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by