Optimizing code which uses repmat and bsxfun

Question

I have written below code to calculate a matrix:

 vec0=repmat(vec,1,9);
 triw = (bsxfun(@times,vecO(1,:)',yc1)-bsxfun(@times,vecO(2,:)',xc1)).*(bsxfun(@times,vecO(2,:)',yc1)+bsxfun(@times,vecO(1,:)',xc1));

vec is a 2-by-900 matrix and xc1 and yc1 are 8100-by-900. I use this code in a loop. It is very slow, so I'd like to improve its performance. How would I do this?

Answer 1

I was able to speed things up by around 30-40% by reorganizing the calculations into 2 calls to bsxfun instead of 4:

triw = bsxfun(@times, prod(vec0).', yc1.^2-xc1.^2)-...
       bsxfun(@times, diff(vec0.^2).', xc1.*yc1);

Note that I also used the array transpose operator .' instead of the complex conjugate transpose operator '. The first simply reorganizes the array without modifying the values, while the second can give you different results if you're dealing with complex data.

Here's the code I used to compare the two approaches:

% Random test data:
vec0 = rand(2, 8100);
xc1 = rand(8100, 900);
yc1 = rand(8100, 900);

% Define anonymous functions to test:
fcn1 = @(v, x, y) (bsxfun(@times, v(1, :).', y)-bsxfun(@times, v(2, :).', x)).*...
                  (bsxfun(@times, v(2, :).', y)+bsxfun(@times, v(1, :).', x));
fcn2 = @(v, x, y) bsxfun(@times, prod(v).', y.^2-x.^2)-...
                  bsxfun(@times, diff(v.^2).', x.*y);

% Test the mathematical accuracy:
triw1 = fcn1(vec0, xc1, yc1);
triw2 = fcn2(vec0, xc1, yc1);
max(abs(triw1(:)-triw2(:)))

ans =
     4.440892098500626e-16

% Time the results:
t1 = timeit(@() fcn1(vec0, xc1, yc1))

t1 =
   0.107775908242267  % seconds

t2 = timeit(@() fcn2(vec0, xc1, yc1))

t2 =
   0.068403928901861  % seconds

The maximum absolute difference between the two results was on the order of the floating-point relative accuracy for double-precision numbers, so effectively no difference.