Julia: `hypot` performance questionable

[KlausC]

I just made massive use of hypot(a, b) with real arguments. I discovered an interesting performance bottleneck.
I am aware of #31922 and #33224, which both focus on accuracy.
I want to present the following measurements:

julia> @benchmark hypot(1.0, 2.0)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     18.065 ns (0.00% GC)
  median time:      18.232 ns (0.00% GC)
  mean time:        18.726 ns (0.00% GC)
  maximum time:     60.926 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     997

julia> @benchmark MatrixAlgebra._hypot(1.0, 2.0)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.830 ns (0.00% GC)
  median time:      1.843 ns (0.00% GC)
  mean time:        1.890 ns (0.00% GC)
  maximum time:     10.976 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

That is a runtime factor of 10 compared to the (maybe naive) implementation:

 function _hypot(x, y)
     s = sqrt(abs2(x) + abs2(y))
     _isclean(s) ? s : hypot(x, y)
 end
 function _isclean(a::T) where T<:AbstractFloat
     isfinite(a) &&
    a >= sqrt(4 * floatmin(T))
end

Comparing the accuracy of both implementations gives:

julia> foo(10^6, Float64, hypot)
(4.005836504709864e-17, 4.7168514470489265e-17, 1.168161969438445e-16)

julia> foo(10^6, Float64, MatrixAlgebra._hypot)
(4.6543949325418543e-17, 5.704669978693873e-17, 2.108974170018499e-16)

julia> function foo(n, T, f)
           s3 = s2 = s1 = zero(T)
           for i = 1:n
               a, b = T.(randn(2))
               h = hypot(big(a), big(b))
               d = T(abs(h - f(a, b)) / h)
               s3 = max(s3, d)
               s2 = hypot(s2, d)
               s1 += d
           end
           s1 / n, s2 / sqrt(n), s3
       end
foo (generic function with 3 methods)

The return values of foo are the mean values of the relative deviations from the accurate values in 1-norm, 2-norm, and Inf-norm.
I want to discuss, if those deviations in accuracy justify the performance burden.

[longemen3000]

I like this fast track, the only thing that comes to my mind Is if the ulp error is enough

[kimikage]

I am concerned about the degradation of accuracy, too.

However, as for the Float32, I think the following specialization may be available (for the intermediate variable s) in most cases.

hypot_f32(a::Float32, b::Float32) = Float32(sqrt(Float64(a)^2 + Float64(b)^2))

(I have not identified the cases which would lead to inaccuracy yet.)

[kimikage]

BTW, the current implementation of hypot uses Base.Math.FMA_NATIVE. (cf. #33011)

[KlausC]

Can anybody run the benchmarks on a machine, with FMA_NATIVE == true? For example as of Intel I5 Generation 4.

[kimikage]

The following are rough (i.e. unfair) benchmarks.

julia> versioninfo()
Julia Version 1.4.2
Commit 44fa15b150* (2020-05-23 18:35 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-8565U CPU @ 1.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, skylake)

julia> @benchmark hypot.(x, y) setup=(x=randn(1000); y=randn(1000))
BenchmarkTools.Trial:
  memory estimate:  7.94 KiB
  allocs estimate:  1
  --------------
  minimum time:     12.600 μs (0.00% GC)
  median time:      13.300 μs (0.00% GC)
  mean time:        14.513 μs (0.81% GC)
  maximum time:     1.190 ms (98.44% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> @benchmark hypot_fma.(x, y) setup=(x=randn(1000); y=randn(1000))
BenchmarkTools.Trial:
  memory estimate:  7.94 KiB
  allocs estimate:  1
  --------------
  minimum time:     8.899 μs (0.00% GC)
  median time:      9.100 μs (0.00% GC)
  mean time:        9.763 μs (1.07% GC)
  maximum time:     1.063 ms (98.65% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> foo(10^6, Float64, hypot_fma)
(4.004969340992502e-17, 4.715230149802108e-17, 1.1094378857772951e-16)

[KlausC]

would also be interested in the times on that machine, using _hypot from my post.

[kimikage]

Please wait a minute. :thinking:

Does this eps fail the constant folding?
https://github.com/JuliaLang/julia/blob/d6d5208d66be637543232bec4b0c79b2811ea1c6/base/math.jl#L662

cf. https://github.com/JuliaLang/julia/pull/31922#issuecomment-489678018
cc: @cfborges, @simonbyrne

[simonbyrne]

It looks like it does...

julia> f(T) = eps(sqrt(floatmin(T)))
f (generic function with 1 method)

julia> @code_llvm f(Float64)

;  @ REPL[9]:1 within `f'
define double @julia_f_17254(%jl_value_t addrspace(10)*) {
top:
; ┌ @ float.jl:745 within `eps'
   %1 = call double @julia_ldexp_918(double 0x3CB0000000000000, i64 -511)
; └
  ret double %1
}

[simonbyrne]

This seems to work correctly:

julia> g(T) = eps(T)*sqrt(floatmin(T))
g (generic function with 1 method)

julia> @code_llvm g(Float64)

;  @ REPL[21]:1 within `g'
define double @julia_g_17323(%jl_value_t addrspace(10)*) {
top:
  ret double 0x1CC0000000000000
}

[KlausC]

PR #36365 improves the relation from 10-fold to 7-fold on my machine (without FMA_NATIVE). That may go down to 5-fold on FMA machines.
Question is: do we want to pay the price for 1 bit of accuracy in hypot?

[cfborges]

I had noticed that evaluating those constants seems to slow things down but did not know why. Glad someone has a notion about that.
It's clear that a code that avoids unnecessary floating exceptions and attempts to correctly round will be slower. My goal with that algorithm was simply to get to correct rounding since that is what the IEEE standard suggests (but does not require) for hypot. That said, I think for most applications where speed is critical it's better to just use h = sqrt(x*x+y*y) since this can never be off by more than one ulp (as opposed to the Julia 1.1 version of hypot which could be off by two). The only problem with this approach is that the user needs to be sure that there are no extreme arguments since those could lead to unecessary overflow/underflow.

[kimikage]

micro optimization...

@inline function hypot(x::T, y::T) where T<:AbstractFloat
    absx, absy = abs(x), abs(y)

    # Return Inf if either or both inputs is Inf (Compliance with IEEE754)
    isinf(absx) && return absx
    isinf(absy) && return absy

    # Order the operands (ax >= ay)
    ax = ifelse(absx < absy, absy, absx)
    ay = ifelse(absx < absy, absx, absy)
    ...

Edit:
I applied the suggestion by @cfborges (https://github.com/JuliaLang/julia/issues/36353#issuecomment-649090655)
They seem to have flipped while I was struggling with NaN. The intention is to use the hardware min/max instructions (without @fastmath). So my bug should have no effect on "speed".

[simonbyrne]

@kimikage can you open a pull request?

[kimikage]

Perhaps I can, but I don't think it's worth the separate PR now. We still need the discussion about #33011, the approach of OP and the Float32 specialization.

[cfborges]

One minor comment I think

# Order the operands
    ax = ifelse(absx < absy, absx, absy)
    ay = ifelse(absx < absy, absy, absx)

would actually need to be

# Order the operands
    ax = ifelse(absx < absy, absy, absx)
    ay = ifelse(absx < absy, absx, absy)

Since the goal is to have x > y.

[simonbyrne]

I'll leave this open since my fix was only partial