Can you do this calculation with bsxfun instead of repmat?

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

Answer 1

[Guillaume]

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

[Henrik]

Cool, I will try it out, thanks for the help!

[Guillaume]

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]

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.