Element-wise raise an array to an integer power

Hi everybody,

I found that the element-wise multiplication of a large array (a.*a.*a) is much faster than the equivalent element-wise power function (a.^3). It seems that the only exception is squaring:

>> a = rand(300,300,300);

>> tic,a.^2;toc

Elapsed time is 0.125585 seconds.

>> tic,a.^3;toc

Elapsed time is 2.304034 seconds.

>> tic,a.^12;toc

Elapsed time is 2.328101 seconds.

while:

>> tic,a.*a;toc

Elapsed time is 0.147440 seconds.

>> tic,a.*a.*a;toc

Elapsed time is 0.214431 seconds.

>> tic,(a.*a.*a).^2.^2;toc

Elapsed time is 0.505067 seconds.

Is there an alternative Matlab element-wise power function which factorizes integer exponents and performs element-wise multiplications instead of exponentiations, whenever preferable?

Thanks a lot

Bernhard

[dpb]

That is disappointing, indeed.

Appears no optimization has been implemented and AFACT there's no alternative builtin in that does.

Looks like "roll your own" time; might check the File Exchange and see if somebody has already done the heavy lifting.

[Walter Roberson]

Pretty sure there is appropriate logic in the file exchange

[dpb]

Another of those bizarre mysteries why TMW wouldn't have done when first wrote the function for such a fundamental operation with the solution so well known. Generally runtime libraries of compilers do; seems as though would get "for free" if just used called it in the bowels.

Answer 1

[Jan]

What about:

• https://www.mathworks.com/matlabcentral/fileexchange/44574-elementwise-power
• https://www.mathworks.com/matlabcentral/fileexchange/44575-matrix-power (matrix multiplications)
• https://www.mathworks.com/matlabcentral/fileexchange/25782-mpower2-a-faster-matrix-power-function (matrix multiplication)

Answer 2

[Walter Roberson]

On my system, with R2018a, timeit shows less discrepancy:

>>  a = rand(300,300,300);
>> timeit(@() a.^2, 0)
ans =
            0.060887501623
>> timeit(@() a.^3, 0)
ans =
            0.544000124623
>> timeit(@() a.^12, 0)
ans =
            0.535817926623
>> timeit(@() (a.*a.*a).^2.^2, 0)
ans =
            0.135812291623

Now, the A.^B operation is defined in terms of exp(log(A)*B) (somewhere... I don't see the reference right this moment.) So I tested max(max(max(abs(exp(log(a)*12) - a.^12)))) -- and it isn't 0, it is eps/2 . So whatever MATLAB is doing internally is not the log implementation, and from the timings is not an integer factorization based implementation either. (I also created a random complex matrix and tested the results of various factorizations, and none of them match bit for bit.)