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?

*No products are associated with this question.*

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.

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.

Opportunities for recent engineering grads.

## 0 Comments