## Generate a vector of powers more quickly using cumprod in MATLAB

While working on someone’s MATLAB code today there came a point when it was necessary to generate a vector of powers. For example, [a a^2 a^3….a^10000] where a=0.999

a=0.9999; y=a.^(1:10000);

This isn’t the only way one could form such a vector and I was curious whether or not an alternative method might be faster. On my current machine we have:

>> tic;y=a.^(1:10000);toc Elapsed time is 0.001526 seconds. >> tic;y=a.^(1:10000);toc Elapsed time is 0.001529 seconds. >> tic;y=a.^(1:10000);toc Elapsed time is 0.001716 seconds.

Let’s look at the last result in the vector y

>> y(end) ans = 0.367861046432970

So, 0.0015-ish seconds to beat.

>> tic;y1=cumprod(ones(1,10000)*a);toc Elapsed time is 0.000075 seconds. >> tic;y1=cumprod(ones(1,10000)*a);toc Elapsed time is 0.000075 seconds. >> tic;y1=cumprod(ones(1,10000)*a);toc Elapsed time is 0.000075 seconds.

soundly beaten! More than a factor of 20 in fact. Let’s check out that last result

>> y1(end) ans = 0.367861046432969

Only a difference in the 15th decimal place–I’m happy with that. What I’m wondering now, however, is will my faster method ever cause me grief?

This is only an academic exercise since this is not exactly a hot spot in the code!

You used a=0.9999, not a=0.999.

If I do this in Euler Math Toolbox, I get a three times better result for a^(…) and a two times worse result for cumprod(…), reducing the factor to about 3. Actually, I do not know, why Matlab gets a different factor. By the way, C is about two times faster than cumprod.

>a=0.9999; n=10000;

>tic; loop 1 to 1000; a^(1:n); end; toc;

Used 0.558 seconds

>tic; loop 1 to 1000; cumprod(a*ones(1,n)); end; toc;

Used 0.118 seconds

>longest a^(1:n)[-1]

0.3678610464329705

>longest cumprod(a*ones(1,n))[-1]

0.3678610464329692

>function tinyc test (a,n) …

$ARG_DOUBLE(a);

$ARG_INT(n);

$RES_DOUBLE_MATRIX(m,1,n);

$int i;

$m[0]=a;

$for (i=1; itic; loop 1 to 1000; test(0.9999,10000); end; toc;

Used 0.077 seconds

>>You used a=0.9999, not a=0.999.

You are quite right. Thank you for pointing out the typo (now fixed) and thanks for trying it out in Euler Math Toolbox. I see that you’ve used the tiny C compiler. Is it as easy to use other compilers in Euler?

In a comment to the Google+ post about this article — https://plus.google.com/108930714356162091073/posts/hu9Du4dBWwY — Matt Tearle said ‘I tried exp(log(a)*(1:10000)) — it was faster than .^ but still a bit slower than the unfortunately-named function.’

You can use any compiler, which produces a DLL.

You can also use Python. Let me try.

>function python test (a,n) …

$ v=[a] * n;

$ for i in range(2,n):

$ v[i]=v[i-1]*a

$ return v

$ endfunction

>tic; loop 1 to 1000; test(0.999,1000); end; toc

Used 0.187 seconds

0.187

This a 2.5 the speed of C. Not bad, and maybe there are better ways in Python.

Could it be that .^ uses the binary method for exponentials? That would explain, why it is slower than exp(n*log(a)).

Using zeros(1,n)+a with cumprod() is ~20% faster than using ones(1,n)*a:

>> tic, for idx=1:1000, clear y; y=cumprod(ones(1,n)*a); end, toc

Elapsed time is 0.053967 seconds.

>> tic, for idx=1:1000, clear y; y=cumprod(zeros(1,n)+a); end, toc

Elapsed time is 0.042769 seconds.

The result sounds very logical. Afterall, cumprod or cumsum requires less number of operations compared to a power (N-1 vs [N*(N+1)-1]/2, where N is the length of the vector).

@Yair Altman This could be even slightly faster, if you use

>> m = []; m(10000) = 0;

instead of

>> m = zeros(1, n)