Computing time for .^(1/2) or .^(1/p)

Hi,

I realised that computing the square root element-wise of a matrix using .^(1/2) is faster than computing the pth root whatever the value of p :

I found (with the provided code below):

Averaged time for A.^(1/2): 0.0337 s

Averaged time for A.^(1/3): 0.0972 s

Averaged time for A.^(1/5): 0.0946 s

Does anyone know why ? (I am clearly not an expert in algorithm complexity and optimization, but I'd like to know if there is a special case of the square root implemented in the function (.^), and eventually to understand how it works in a simple way).

Thanks in advance :) !

clear all
close all
clc

% Data
A = randn(1024,1024);


% Averaged time to Compute the square root of each element in A
t_squareroot = 0;
for u = 1 : 50
    tic
    B = A.^(1/2);
    t_squareroot = t_squareroot + toc/50;
end
disp('Averaged time for A.^(1/2):')
disp(t_squareroot)


% Averaged time to Compute the cubic root of each element in A
t_cubicroot = 0;
for u = 1 : 50
    tic
    C = A.^(1/3);
    t_cubicroot = t_cubicroot + toc/50;
end
disp('Averaged time for A.^(1/3)')
disp(t_cubicroot)


% Averaged time to Compute the 5th root of each element in A
t_5throot = 0;
for u = 1 : 50
    tic
    D = A.^(1/5);
    t_5throot = t_5throot + toc/50;
end
disp('Averaged time for A.^(1/5)')
disp(t_5throot)

Answer 1

[Stephane Dauvillier]

Explanation in french here

In order to compute the square root of any number you only need to be able to compute the square root of the mantisse (which is between 1 and 2). You only have to subtract 1 to the exposant which is absolutly costless.

There is no such thing for the other power.

[James Tursa]

I don't see how the "subtract 1" part works for the exponent ... doesn't seem like it would do the trick. Can you explain and/or provide an example?

[Stephane Dauvillier]

Sorry, not subtract by one but shift one time (meaning divide by 2)

[James Tursa]

And what to do if the exponent is not a multiple of 2?

The point I am driving at is it is not quite as simple as you originally stated. You first need a method that can do a quick square root of the mantissa (faster than a generic method for an arbitrary size number) or the whole exercise is pointless. Then you need to account for exponents that are not divisible by 2.

You could do something similar for other roots but I don't know if anyone has ever bothered. E.g., you could develop a rational function approximation to the cube root over the mantissa range, and then some bit manipulation to fix up the exponent, etc.

[Walter Roberson]

Subtract 1 does not sound right, I agree.

https://pdfs.semanticscholar.org/5060/4e9aff0e37089c4ab9a376c3f35761ffe28b.pdf

https://en.wikipedia.org/wiki/Methods_of_computing_square_roots

[Stephane Dauvillier]

If exponent is not divided by 2 then

NewEponent =(oldExponent-1)/2;

NewMantissa=sqrt(2*oldMantissa)

[Maxime Polichetti]

Thank you for answering :). I can understand now how square root could be faster than others.

One last thing, I can't find in the matlab documentation, that they use such a "trick" for it (I tried F1 for sqrt, power, nthroot, functions and their corresponding m version)... but nothing in the "More About" Section. Is it a "confidential" information ?

Thanks again :) !

[Walter Roberson]

In situations such as these, it is not necessarily the case that similar speed-ups could not be found for other specific roots, but rather that they may need custom coding for each different root, and so there become questions of whether it is worth writing any particular one -- is the programming/debugging time and code maintenance time worth the extra performance compared to the amount of time that particular exponent is used?

sqrt() gets used a lot, so it is a natural. Cube root... maybe . 4th root... I do see that used, but I probably see random non-integer roots used more often.

Answer 2

[James Tursa]

I am unaware of any documentation on how TMW has chosen to implement their rooting/power algorithms, but it would not surprise me that they use a custom fast algorithm for sqrt vs a slower more generic algorithm for other powers & roots. E.g., one could use a rational function approximation for the mantissa combined with direct exponent bit manipulation for square root, and then something else for all the other powers. (This is just an example of what could be done ... not necessarily my guess at what the library code that MATLAB uses does do).

Answer 3

[Walter Roberson]

https://www.felixcloutier.com/x86/sqrtpd (hardware square root)

[KALYAN ACHARJYA]

oK sir