isapprox: test max(atol,rtol*...) rather than atol+rtol*...

[stevengj]

Ever since the isapprox function was merged in #3408, we have been testing something like norm(x-y) <= atol + rtol*max(norm(x), norm(y)). This PR changes that to norm(x-y) <= max(atol, rtol*max(norm(x), norm(y))).

Rationale: by adding the two tolerances when atol ≠ 0, we were effectively increasing the tolerance for numbers with atol close to rtol*norm(x), which doesn't seem desirable (though it is not terrible). Python 3's new is_close function introduced in PEP 485 uses the max behavior, which seems a bit more sensible.

As a practical matter, I seriously doubt this will break any real code, but technically it is a breaking change.

(As an aside, the PEP contains this remark that should be gratifying to Julians: In theory, it should work for any type that supports abs() , multiplication, comparisons, and subtraction. However, the implementation in the math module is written in C, and thus can not (easily) use python's duck typing.)

[stevengj]

The only argument I can think of for the old behavior is performance: + is cheaper than max. In a quick test of x = 1.0; y = 1.001; @btime isapprox($x, $y);, I find 9.05ns in 0.6 vs 9.35ns for the new method in 0.7, a 3% slowdown. I seriously doubt this matters in any real application, where isapprox is rarely a critical inner-loop function.

[stevengj]

Hmm, it seems I did the wrong benchmark comparisons by comparing to 0.6. In 0.7 master (just before this PR), I get 7.85ns for that benchmark, a big improvement over 0.6. Compared to that, this PR is a 20% slowdown, which is surprising.

(Anyone know why isapprox got so much faster in 0.7 than in 0.6?)

I still tend to think it is worth it, but the performance rationale for the old behavior is stronger than I expected.

[StefanKarpinski]

No idea about the performance change, but I agree that it seems unlikely that isapprox is performance-critical anywhere. It tends to be used to verify values at the end of a computation. This change does seem to be more correct, so it seems like a good change to make now.

[stevengj]

Travis failure is "The job exceeded the maximum time limit for jobs, and has been terminated.", so presumably unrelated.

[stevengj]

Adding this to 1.0 to make sure it is decided by then.

[stevengj]

Rebased.

[StefanKarpinski]

I'm all for it. @jiahao, @andreasnoack? Trying to think of others who haven't already thumbed this up who might care about this issue...

[andreasnoack]

I think this is an improvement but it seems that this way of testing is basically broken when we allow both tolerances to be non-zero. E.g.

julia> x, y = 1e8, 1e8 + 0.5;

julia> atol = 0.1;

julia> abs(x-y) <= max(atol, Base.rtoldefault(x, y)*max(abs(x), max(y)))
true

julia> abs(x-y) <= atol
false

Maybe it would be a better API to have tol and relative::Bool keywords. In any case, this PR is better than what we have now so these considerations shouldn't block it.

[StefanKarpinski]

If having both be non-zero is already broken it would seem better to get the design right rather than change this and then also deprecate it. Although I suppose that changing this could pre-emptively find some cases where an API change would break overly loose tests anyway.

[Sacha0]

Thoughts from a brief chat with @andreasnoack. (Thanks Andreas! :) )

Three designs seem reasonable depending on what is useful in practice:

A) If simultaneously enforcing absolute and relative tolerances isn't useful in practice, then some form of @andreasnoack's suggestion seems wise: Allow the user to specify either an absolute or relative tolerance, but not both.

B) If simultaneously enforcing absolute and relative tolerances is useful in practice, then the following design seems reasonable:

1. If the user specifies only an absolute tolerance, guarantee only that absolute tolerance
2. If the user specifies only a relative tolerance, guarantee only that relative tolerance
3. If the user specifies both absolute and relative tolerances, guarantee both tolerances independently

... i.e. min seems like a reasonable way to combine tolerances when the user specifies both.

C) If both simultaneously enforcing absolute and relative tolerances and flexibility in their combination are useful in practice, then the following design seems reasonable:

1. If the user specifies only an absolute tolerance, guarantee only that absolute tolerance
2. If the user specifies only a relative tolerance, guarantee only that relative tolerance
3. If the user specifies both absolute and relative tolerances, require an additional keyword (e.g. combine) specifying how the tolerances should be combined (a function, e.g. min, max, or +)

Chances are scheme (C) is unnecessarily complex. One way to mitigate the complexity would be to make combine default to e.g. min, yielding the behavior of scheme (B) by default while providing flexibility if/when need be. Best!

[stevengj]

You occasionally need both tolerances. Typically in cases where you are evaluating a function at lots of points and what you mostly want is a relative tolerance, but you don't want it to unexpectedly break at points where the result happens to be nearly zero.

Note that this implies max, not min. The min approach seems inherently flawed to me.

[stevengj]

That is, I see Andreas' "flaw" as a feature. You want it to return true if either test is satisfied, not necessarily both, in any use case I'm aware of.

[stevengj]

Another use case of this sort comes up in convergence tests for quadrature, where you want to halt if the estimated integral isapprox equal to a lower-order estimate. You usually want a relative tolerance, but it can be useful to specify an absolute tolerance too just in case the integral is nearly zero, to make sure it terminates. You want it to terminate when either tolerance is satisfied. (This is typical of lots of iterative methods where you specify multiple convergence criteria.)

[andreasnoack]

You want it to return true if either test is satisfied

That makes sense but I think it is a problem that approx(x, y, atol = something) suggests that you test the absolute difference and it is not obvious that a hidden default value for the relative tolerance might make the test pass. That could be fixed by setting the default value for rtol=0 and only set it to rtoldefault(x,y) when atol==rtol==0.

[StefanKarpinski]

So that suggests C with combine=max?

[Sacha0]

@stevengj, thank you for the super-useful concrete use cases! :)

So that suggests C with combine=max?

Yes, perhaps then scheme (C) is in order :). isapprox(x, y, atol = foo) would work as @andreasnoack desires, and isapprox(x, y, atol = foo, rtol = bar, combine = max) (with whichever default well documented) would cover @stevengj's use case? Best!

[stevengj]

I'm skeptical of the need for a combine argument, but @andreasnoack's suggestion of a default rtol=atol=0 (setting rtoldefault if both are zero) makes sense to me.

[stevengj]

We could also add a combine argument later if the need arises, since it wouldn't be breaking as long as the default is max. But honestly, I've never heard of a need for this.

[Sacha0]

Cheers! Sounds like scheme (B) with max rather than min, leaving the option to extend to scheme (C) if/when use cases appear? That is,

1. If the user does not specify any tolerances, guarantee only some default relative tolerance
2. If the user specifies only an absolute tolerance, guarantee only that absolute tolerance
3. If the user specifies only a relative tolerance, guarantee only that relative tolerance
4. If the user specifies both absolute and relative tolerances, guarantee the maximum of those tolerances

? Best!

[stevengj]

Updated to the proposed semantics (default rtol=0 if atol>0 is supplied).

[tkelman]

enough packages are calling the two-argument version (despite it not being exported) that it should probably be deprecated

[stevengj]

Okay, I just pushed a deprecation.

[stevengj]

Looks like an unrelated Travis failure (read end of file in the distributed.jl test).

[amitmurthy]

The distributed.jl failure both here

	From worker 2:	signal (11): Segmentation fault
	From worker 2:	while loading no file, in expression starting on line 0
	From worker 2:	jl_typemap_assoc_by_type_ at /home/travis/build/JuliaLang/julia/src/typemap.c:571
	From worker 2:	jl_mt_assoc_by_type at /home/travis/build/JuliaLang/julia/src/gf.c:1041
	From worker 2:	jl_lookup_generic_ at /home/travis/build/JuliaLang/julia/src/gf.c:1897 [inlined]
	From worker 2:	jl_apply_generic at /home/travis/build/JuliaLang/julia/src/gf.c:1926
	From worker 2:	jl_apply at /home/travis/build/JuliaLang/julia/src/julia.h:1448 [inlined]
	From worker 2:	start_task at /home/travis/build/JuliaLang/julia/src/task.c:268
	From worker 2:	Allocations: 1321563 (Pool: 1321465; Big: 98); GC: 7

and in #22998

signal (11): Segmentation fault: 11
while loading /private/tmp/julia/share/julia/test/distributed_exec.jl, in expression starting on line 1333
jl_lookup_generic_ at /Users/travis/build/JuliaLang/julia/src/gf.c:1873
jl_apply_generic at /Users/travis/build/JuliaLang/julia/src/gf.c:1926
jl_apply at /Users/travis/build/JuliaLang/julia/src/./julia.h:1448 [inlined]
start_task at /Users/travis/build/JuliaLang/julia/src/task.c:268
Allocations: 27875061 (Pool: 27872524; Big: 2537); GC: 76

are a segfault in gf.c .

Has it been noticed in any other PRs?

Maybe @vtjnash / @yuyichao can help?

[JeffBezanson]

Bump from triage. Needs a rebase then good to go.

[stevengj]

Rebased.

[Jutho]

Nice to discover why the tests I wrote earlier today suddenly fail. I was falling in the class of problems that @stevengj described, where I test a lot of non-zero numbers, but occasionally the result is zero so I need an atol on top of a rtol. I think this is indeed the most common use case for ≈.