Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

To resolve issues starting MATLAB on Mac OS X 10.10 (Yosemite) visit: http://www.mathworks.com/matlabcentral/answers/159016

Is there a better vectorization technique?

Asked by Kalamaya on 29 Aug 2012

I am trying to see if there are other ways of coding this code sample more efficiently. Here, y is an 1xM matrix, (say, 1x1000), and z is an NxM matrix, (say, 5x1000).

   mean(ones(N,1)*y.^3 .* z,2)

This code works fine, but I worry of N increases a lot, that the:

   ones(N,1)*y.^3 

might get too wasteful and make everything slow down.

Thoughts?

0 Comments

Kalamaya

Products

No products are associated with this question.

2 Answers

Answer by John D'Errico on 29 Aug 2012
Accepted answer

Its not THAT terrible for a matrix that small. However, you can gain from the use of bsxfun in many instances. Here, the matrices are simply too small to really gain anything.

>> N = 5;M =1000;
>> y = rand(1,M);
>> z = rand(N,M);
>> mean(ones(N,1)*y.^3 .* z,2)
ans =
      0.12412
      0.11669
      0.12102
      0.11976
      0.12196
>> mean(bsxfun(@times,y.^3,z),2)
ans =
      0.12412
      0.11669
      0.12102
      0.11976
      0.12196
>> z*y.'.^3/M
ans =
      0.12412
      0.11669
      0.12102
      0.11976
      0.12196

As you can see, all three solutions return the same result. All are equally valid.

Now I'll compare the times required.

>> timeit(@() mean(ones(N,1)*y.^3 .* z,2))
ans =
   0.00023018
>> timeit(@() mean(bsxfun(@times,y.^3,z),2))
ans =
   0.00026829
>> timeit(@() z*y.'.^3/M)
ans =
   0.00016594

As I said, you don't gain much. In fact, bsxfun does not gain at all, and is a bit slower. But you can gain a bit, if you re-write the expression into the third form I've posed. Not much, but a bit.

2 Comments

Kalamaya on 29 Aug 2012

I thought about using bsxfun. How is the timing if the N was very large, say, 2000 or something like that? (timeit is your own function?)

John D'Errico
Answer by Jan Simon on 29 Aug 2012
Edited by Jan Simon on 30 Aug 2012

The power() operator is more expensive than the matrix multiplication. Therefore an explicit multiplication saves time:

M = 1000;
N = 5;
y = rand(1, M);
z = rand(N, M);
tic; for i=1:100, a = mean(ones(N,1) * y .^ 3 .* z, 2); end; toc
% Elapsed time is 0.036668 seconds.
tic; for i=1:100, a = z * y.' .^ 3 / M; end; toc
% Elapsed time is 0.026818 seconds.
tic; for i=1:100, a = z * (y .* y .* y)' / M; end; toc
% Elapsed time is 0.001327 seconds.

[EDITED] If the resulting array is very large, a multiplication is faster than a division, but the result can differ due to rounding:

a = z * (y .* y .* y)' * (1 / M);

For the small [5x1] array in the example this does not matter.

0 Comments

Jan Simon

Contact us