Can you do this calculation with bsxfun instead of repmat?

1 view (last 30 days)
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

Accepted Answer

Guillaume
Guillaume on 26 Sep 2014
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 Comments
Guillaume
Guillaume on 27 Sep 2014
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 on 27 Sep 2014
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!