fixed hypot(x::Number...) for under/overflow #27251

johnfgibson · 2018-05-24T19:51:29Z

For issue #27141. The code is taken from vecnorm functions in LinearAlgebra/generic.gl with a few replacements of norm with abs.

(py36) gibson@sophist$ ./julia 
               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: https://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.7.0-DEV.5214 (2018-05-24 11:25 UTC)
 _/ |\__'_|_|_|\__'_|  |  Commit 14dbdf509a* (0 days old master)
|__/                   |  x86_64-suse-linux

julia> hypot(1e-300, 0.0, 0.0)
1.0e-300

julia> hypot(1e-300 + 0im, 0.0, 0.0)
1.0e-300

stevengj · 2018-05-25T07:50:05Z

base/math.jl

+        vnorm = abs(v)
+        maxabs = ifelse(isnan(maxabs) | (maxabs > vnorm), maxabs, vnorm)
+    end
+    maxabs = float(maxabs)


Type instability?

Yes. I'll see if I can fix that by determining the least common supertype of the elements of the Tuple x and promoting to that over the iterations.

Why not using another variable?

Using another variable for line 517 is good, but @code_warntype showed more type instability due to the variability of types within the Tuple argument x, which xp = promote(x...) fixes. I'm just about ready to commit a type-stable version.

StefanKarpinski · 2018-05-25T17:35:47Z

base/math.jl

-    # compute infnorm x (modeled on generic_vecnormMinusInf(x) in LinearAlgebra/generic.gl)
-    (v, s) = iterate(x)::Tuple
+
+    xp = promote(x...)


Perhaps better to just make the signature hypot(x::T...) where T<:Number and then have a separate method for hypot(x::Number...) = hypot(promote(x...))? Otherwise the compiler has to prove that this is a no-op to eliminate it—which it probably can—but two methods seems clearer anyway.

Sounds good. I'll commit a slight revision of that: hypot(x::Number...) = hypot(promote(x...)...) correctly passes off to function hypot(x::T...) where T<:Number,

whereas hypot(x::Number...) = hypot(promote(x...)) results in

julia> hypot(1.0, 0, 0) ERROR: MethodError: no method matching hypot(::Tuple{Float64,Float64,Float64})

Thanks for all the help. I'm learning a lot.

Yes, that's what I meant :)

With the above revisions make test fails with the following. I see from NEWS.md that this has to do with "possibly throw[ing] UndefVarError when accessing static parameters," but I've got no idea how to address it.

Error in testset ambiguous: Test Failed at /home/gibson/gitworking/julia/test/ambiguous.jl:302 Expression: need_to_handle_undef_sparam == Set() Evaluated: Set(Method[hypot(x::T...) where T<:Number in Base.Math at math.jl:510]) == Set(Any[]) ERROR: LoadError: Test run finished with errors in expression starting at /home/gibson/gitworking/julia/test/runtests.jl:59

It seems like the ambiguous case is hypot(). What should the hypotenuse of zero numbers be?

nothing? 0.0? Throwing an error seems most sensible to me. That's what norm([]) does.

In that case the fix would be to have a method that throws an error for zero arguments similar to what norm does.

stevengj · 2018-05-29T19:48:49Z

base/math.jl

+        vnorm = abs(v)
+        maxabs = ifelse(isnan(maxabs) | (maxabs > vnorm), maxabs, vnorm)
+    end
+    maxabsf = float(maxabs)


What is the advantage of this over:

maxabsf = float(maximum(abs, x)) isnan(maxabsf) && return maxabsf

stevengj · 2018-05-29T19:51:09Z

base/math.jl

+            (v, s) = y
+            sum += abs2(v)
+        end
+        return convert(Tfloat, sqrt(sum))


Why write out the loop, instead of just doing

return oftype(sqrt(sum(x -> oftype(maxabsf, abs2(x)), x)), maxabsf)

stevengj · 2018-05-29T19:53:48Z

base/math.jl

+            sum += (abs(v)/maxabsf)^2
+        end
+        return convert(Tfloat, maxabsf*sqrt(sum))
+    end


Similarly, why not just:

return oftype(sqrt(sum(x -> (abs(x)/maxabsf)^2, x)), maxabsf)

I copied generic_vecnorm2 from generic.jl in LinearAlgebra, inlining the call to vecnormInf on the way. I take it your compact code is as efficient. Let me pull it apart to understand it better and verify with some benchmarking.

julia/stdlib/LinearAlgebra/src/generic.jl

Lines 298 to 322 in 1d3de4c

function generic_vecnorm2(x)

maxabs = vecnormInf(x)

(maxabs == 0 || isinf(maxabs)) && return maxabs

(v, s) = iterate(x)::Tuple

T = typeof(maxabs)

if isfinite(_length(x)*maxabs*maxabs) && maxabs*maxabs != 0 # Scaling not necessary

sum::promote_type(Float64, T) = norm_sqr(v)

while true

y = iterate(x, s)

y === nothing && break

(v, s) = y

sum += norm_sqr(v)

end

return convert(T, sqrt(sum))

else

sum = abs2(norm(v)/maxabs)

while true

y = iterate(x, s)

y === nothing && break

(v, s) = y

sum += (norm(v)/maxabs)^2

end

return convert(T, maxabs*sqrt(sum))

end

end

@

Unrolling the loops as in the PR is faster and smaller in memory than calling sum(x -> oftype(maxabsf, abs2(x)), x) etc. For four arguments it's ~~~three orders of magnitude~~~ (oops) thirty times faster. As n gets into the thousands the ratio asymptotes to about two.

julia> x = randn(4); julia> @benchmark hypot_unrolled(x...) BenchmarkTools.Trial: memory estimate: 80 bytes allocs estimate: 5 -------------- minimum time: 127.528 ns (0.00% GC) median time: 128.967 ns (0.00% GC) mean time: 142.552 ns (8.53% GC) maximum time: 55.249 μs (99.72% GC) -------------- samples: 10000 evals/sample: 886 julia> @benchmark hypot_callsum(x...) BenchmarkTools.Trial: memory estimate: 560 bytes allocs estimate: 21 -------------- minimum time: 3.965 μs (0.00% GC) median time: 4.032 μs (0.00% GC) mean time: 4.902 μs (16.56% GC) maximum time: 6.260 ms (99.86% GC) -------------- samples: 10000 evals/sample: 8 julia> versioninfo() Julia Version 0.7.0-DEV.5219 Commit fef8b1fffa* (2018-05-25 17:26 UTC) Platform Info: OS: Linux (x86_64-suse-linux) CPU: Intel(R) Core(TM) i7-3960X CPU @ 3.30GHz WORD_SIZE: 64 LIBM: libopenlibm LLVM: libLLVM-6.0.0 (ORCJIT, sandybridge)

for code

hypot_unrolled(x::Number...) = hypot_unrolled(promote(x...)...) function hypot_unrolled(x::T...) where T<:Number # compute infnorm x (modeled on generic_vecnormMinusInf(x) in LinearAlgebra/generic.gl) (v, s) = iterate(x)::Tuple maxabs = abs(v) while true y = iterate(x, s) y === nothing && break (v, s) = y vnorm = abs(v) maxabs = ifelse(isnan(maxabs) | (maxabs > vnorm), maxabs, vnorm) end maxabsf = float(maxabs) # compute vecnorm2(x) (modeled on generic_vecnorm2(x) in LinearAlgebra/generic.gl) (maxabsf == 0 || isinf(maxabsf)) && return maxabsf (v, s) = iterate(x)::Tuple Tfloat = typeof(maxabsf) if isfinite(length(x)*maxabsf*maxabsf) && maxabsf*maxabsf != 0 # Scaling not necessary sum::promote_type(Float64, Tfloat) = abs2(v) while true y = iterate(x, s) y === nothing && break (v, s) = y sum += abs2(v) end return convert(Tfloat, sqrt(sum)) else sum = (abs(v)/maxabsf)^2 while true y = iterate(x, s) y === nothing && break (v, s) = y sum += (abs(v)/maxabsf)^2 end return convert(Tfloat, maxabsf*sqrt(sum)) end end hypot_callsum(x::Number...) = hypot_callsum(promote(x...)...) function hypot_callsum(x::T...) where T<:Number maxabsf = float(maximum(abs, x)) isnan(maxabsf) && return maxabsf (maxabsf == 0 || isinf(maxabsf)) && return maxabsf if isfinite(length(x)*maxabsf*maxabsf) && maxabsf*maxabsf != 0 # Scaling not necessary return oftype(maxabsf, sqrt(sum(x -> oftype(maxabsf, abs2(x)), x))) else return oftype(maxabsf, maxabsf*sqrt(sum(x -> (abs(x)/maxabsf)^2, x))) end end

Don't benchmark calls involving a splatted global. Try

tst_hypot(x) = hypot(x...) @btime tst_hypot($x)

etcetera

Still a factor of sixteen (I misplaced a decimal when I reported 10^3).

julia> x = randn(4); julia> tst_hypot_unrolled(x) = hypot_unrolled(x...); julia> tst_hypot_callsum(x) = hypot_callsum(x...); julia> @btime tst_hypot_unrolled($x) 126.336 ns (5 allocations: 80 bytes) 1.125466836605211 julia> @btime tst_hypot_callsum($x) 2.070 μs (21 allocations: 560 bytes) 1.125466836605211

hypot_callsum is better in terms of code size and maintainability, but it seems to me that hypot_unrolled (which matches the current PR for base/math.jl's hypot) is more in tune with Julia's aim of optimizing speed and accuracy simultaneously.

An alternative, if no one likes adding this function to base, is to remove it and require LinearAlgebra for variable-argument hypot calls. Then this function could be implemented as a call to vecnorm2, as it was prior to #27141.

stevengj · 2018-06-06T21:10:13Z

I'm not convinced that this function is worth the trouble — is there any other language or common library where hypot is defined for more than two arguments? I'm having trouble thinking of a use case that would not be better off calling norm for a StaticArray or similar. Does any package actually use this method?

giordano · 2018-06-06T22:12:54Z

is there any other language or common library where hypot is defined for more than two arguments?

GSL has a three-argument function, but not a variadic function

StefanKarpinski · 2018-06-06T22:31:07Z

What's the harm in having this since we have a correct, fast implementation. Even having it be a generated function is not so bad: if no one ever calls it with k arguments, then who cares?

johnfgibson · 2018-06-07T02:53:48Z

I implemented this function only because issue #27141 shows the existing version in Base is prone to underflow and overflow, and the suggestion there that adapting the vecnorm2 algorithm from LinearAlgebra would fix it. I'm happy for those with better perspective on the whole language to judge whether this PR is a helpful addition or code bloat in Base. But if the latter, the variadic hypot should be removed. Better not to have the function at all than have it give bad results.

johnfgibson · 2018-06-14T18:16:56Z

I'll close this one out. If I can figure out how deprecation works, I'll submit a new PR that deprecates the variadic hypot.

fixed hypot(x::Number...) for under/overflow

6ccd313

ararslan added maths Mathematical functions needs tests Unit tests are required for this change labels May 24, 2018

johnfgibson added 3 commits May 24, 2018 16:30

Add under/overflow tests for hypot

1af76d9

removed empty testset stub

242c94c

correct underflow test

6959b71

StefanKarpinski requested a review from stevengj May 24, 2018 20:50

stevengj reviewed May 25, 2018

View reviewed changes

now type stable

fef8b1f

StefanKarpinski reviewed May 25, 2018

View reviewed changes

throw error for no arguments; now passes tests

74a0a4c

stevengj reviewed May 29, 2018

View reviewed changes

johnfgibson closed this Jun 14, 2018

stevengj mentioned this pull request Dec 8, 2018

fix hypot with more than two arguments #30301

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

fixed hypot(x::Number...) for under/overflow #27251

fixed hypot(x::Number...) for under/overflow #27251

johnfgibson commented May 24, 2018

stevengj May 25, 2018

johnfgibson May 25, 2018

giordano May 25, 2018

johnfgibson May 25, 2018

StefanKarpinski May 25, 2018 •

edited

Loading

johnfgibson May 25, 2018

StefanKarpinski May 25, 2018

johnfgibson May 25, 2018

StefanKarpinski May 25, 2018

johnfgibson May 25, 2018

StefanKarpinski May 26, 2018

stevengj May 29, 2018

stevengj May 29, 2018 •

edited

Loading

stevengj May 29, 2018

johnfgibson May 29, 2018 •

edited

Loading

johnfgibson May 30, 2018 •

edited

Loading

stevengj May 30, 2018

johnfgibson May 30, 2018 •

edited

Loading

johnfgibson Jun 6, 2018

stevengj commented Jun 6, 2018 •

edited

Loading

giordano commented Jun 6, 2018

StefanKarpinski commented Jun 6, 2018

johnfgibson commented Jun 7, 2018

johnfgibson commented Jun 14, 2018 •

edited

Loading

	function generic_vecnorm2(x)
	maxabs = vecnormInf(x)
	(maxabs == 0 \|\| isinf(maxabs)) && return maxabs
	(v, s) = iterate(x)::Tuple
	T = typeof(maxabs)
	if isfinite(_length(x)maxabsmaxabs) && maxabs*maxabs != 0 # Scaling not necessary
	sum::promote_type(Float64, T) = norm_sqr(v)
	while true
	y = iterate(x, s)
	y === nothing && break
	(v, s) = y
	sum += norm_sqr(v)
	end
	return convert(T, sqrt(sum))
	else
	sum = abs2(norm(v)/maxabs)
	while true
	y = iterate(x, s)
	y === nothing && break
	(v, s) = y
	sum += (norm(v)/maxabs)^2
	end
	return convert(T, maxabs*sqrt(sum))
	end
	end

fixed hypot(x::Number...) for under/overflow #27251

fixed hypot(x::Number...) for under/overflow #27251

Conversation

johnfgibson commented May 24, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

StefanKarpinski May 25, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj May 29, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

johnfgibson May 29, 2018 • edited Loading

Choose a reason for hiding this comment

johnfgibson May 30, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

johnfgibson May 30, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj commented Jun 6, 2018 • edited Loading

giordano commented Jun 6, 2018

StefanKarpinski commented Jun 6, 2018

johnfgibson commented Jun 7, 2018

johnfgibson commented Jun 14, 2018 • edited Loading

StefanKarpinski May 25, 2018 •

edited

Loading

stevengj May 29, 2018 •

edited

Loading

johnfgibson May 29, 2018 •

edited

Loading

johnfgibson May 30, 2018 •

edited

Loading

johnfgibson May 30, 2018 •

edited

Loading

stevengj commented Jun 6, 2018 •

edited

Loading

johnfgibson commented Jun 14, 2018 •

edited

Loading