[julia-users] Re: sqrt(abs2(x)) Vs abs(x)

[Steven G. Johnson]

On Monday, May 18, 2015 at 10:51:48 AM UTC-4, antony schutz wrote:
>
> Hello, 
>
> I just discover abs2() function :v 
>
> And I have to say that i'm really surprised that: 
> sqrt(abs2( x ))  is quicker than  abs( x ) 
>

Note that the abs(x) function is more robust because it avoids spurious 
overflow when x is large.

For example, if x=1e300, then abs(x) correctly returns x, but abs2(x) is 
Inf (infinity, i.e. it overflowed the largest representable 
double-precision floating-point number), so sqrt(abs2(x)) = sqrt(Inf) = Inf.

For benchmarking, I used:

       function foo!(Y,X,n)
          for i = 1:length(X), j = 1:n
             @inbounds Y[i] = abs(X[i])
          end
          return Y
       end

       function bar!(Y,X,n)
          for i = 1:length(X), j=1:n
             @inbounds Y[i] = sqrt(abs2(X[i]))
          end
          return Y
       end

For real data, I got:

X = rand(10^3); Y = similar(X);
foo!(Y, X, 10000); bar!(Y, X, 10000); # call once to eliminate JIT cost
@time foo!(Y, X, 10000);
@time bar!(Y, X, 10000);

elapsed time: 0.003777997 seconds (80 bytes allocated)
elapsed time: 0.041755933 seconds (80 bytes allocated)

i.e. abs was much faster.  For complex data, X = rand(Complex128, 10^3), Y 
= Array(Float64, 10^3), I got:

elapsed time: 0.09278846 seconds (80 bytes allocated)
elapsed time: 0.041389585 seconds (80 bytes allocated)

so in this case sqrt(abs2) was faster (though less robust): the abs 
function is costly because it requires both a branch and a rescaling 
(division) in addition to the square root and mul/adds.  See also the 
hypot(x,y) function.