Generate a vector of powers more quickly using cumprod in MATLAB

August 5th, 2013 | Categories: Making MATLAB faster, math software, matlab, programming | Tags:

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!

  1. August 5th, 2013 at 20:42
    Reply | Quote | #1

    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

  2. Mike Croucher
    August 6th, 2013 at 09:20
    Reply | Quote | #2

    >>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?

  3. Mike Croucher
    August 6th, 2013 at 09:23
    Reply | Quote | #3

    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.’

  4. August 6th, 2013 at 13:14
    Reply | Quote | #4

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

  5. August 6th, 2013 at 18:39
    Reply | Quote | #5

    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.

  6. Chandrakanth Terupally
    August 16th, 2013 at 05:14
    Reply | Quote | #6

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

  7. February 19th, 2014 at 08:52
    Reply | Quote | #7

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

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

    instead of

    >> m = zeros(1, n)