MATLAB: Can you do this calculation with bsxfun instead of repmat

bsxfunrepmatspeed

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

Best Answer

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);
Related Question