improve performance of isapprox(::AbstractArray, ::AbstractArray) when rtol = 0

[thchr]

The default value of rtol for isapprox depends on whether atol is zero or not; if atol is nonzero, rtol is zero (#22742).

The "normal-path" isapprox(x, y; atol, rtol check is just norm(x-y) <= max(atol, rtol*max(norm(x), norm(y))). If rtol = 0, however, we are then redundantly computing the norm of both x and y, which can be expensive for arrays (unlike numbers).

This change simply puts a branch in to check iszero(rtol): if so, it avoids computing the norms redundantly. This speeds up isapprox(X, Y; atol ≠ 0) considerably:

julia> A = rand(3, 3)
julia> @btime isapprox($A, $A; atol=1e-11);
  124.592 ns (1 allocation: 128 bytes) # master
julia> @btime isapprox($A, $A; atol=1e-11);
  68.821 ns (1 allocation: 128 bytes) # PR

julia> B = @SMatrix rand(3, 3)
julia> @btime isapprox($B, $B; atol=1e-11);
  14.544 ns (0 allocations: 0 bytes) # master
julia> @btime isapprox($B, $B; atol=1e-11);
  6.900 ns (0 allocations: 0 bytes) # PR

The cost of the extra branch is not nothing, of course, but should be negligible in the context of other costs for most array sizes.

---

There's also a small change to the atol default value: I noticed that care was being taken to initialize rtol to a type consistent with x and y; I figured the same might as well be done for atol.

[LilithHafner]

Nice!

The cost of the extra branch is not nothing, of course, but should be negligible in the context of other costs for most array sizes.

That cost seems negligible for all array sizes to me. That part if this PR looks ready to merge to me.

For the part about atol's type, is there a case where atol=0 is sub-optimal?

[thchr]

After making this PR, I realized that the change to atol's default is not actually unrelated or optional - and probably also not sufficient.
Consider that since there is now a branch and not just a max call, the types of atol and rtol*max(norm(x), ...) need to be the same, otherwise tol will be type-unstable. To exemplify, consider:

julia> f(;x, y) = x>0 ? x : y
julia> @code_warntype f(;x=0, y=0.0)
MethodInstance for (::var"#f##kw")(::NamedTuple{(:x, :y), Tuple{Int64, Float64}}, ::typeof(f))
[...]
Body::Union{Float64, Int64}

So the change to atol's default is necessary. In fact, it is still not precise enough: the type of atol shouldn't just match with rtol's type (it will after this PR), it should match the type of rtol*max(norm(x), ...). This will not be a problem for float eltypes, but it seems to be an issue for e.g. <:Integer eltypes since norm(::AbstractArray{<:Integer} will still return a floating point value (so rtol*max(norm(x), ...) is not <:Integer, but atol's default will be).

This also poses a problem if a call now explicitly has atol=0 - there will then be dynamic dispatch :(.

The problem is complicated by the fact that it is possible to change the used norm via the norm keyword and also by rtol's default already depending on atol. So, it seems we want some way to promote both rtol and atol to a shared type of rtol*max(norm(x), ...) without actually computing norm(x).

Here's a suggestion for that - would love to hear thoughts:

[...]
    d = norm(x - y)
    if isfinite(d)
        T = promote_type(typeof(atol), typeof(rtol*d))
        rtol = convert(T, rtol)
        atol = convert(T, atol)
        tol = iszero(rtol) ? atol : max(atol, rtol*max(norm(x), norm(y))) # <--- now `tol` should be type-stable
        return d <= tol
    else

[thchr]

Okay, that should solve my worry about dynamic dispatch, thanks @LilithHafner.

I think it's fine to keep the more refined atol default though, both for style consistency with efforts for rtol and to make it likely to do same-type <= comparisons.

[LilithHafner]

I'd still like to see an example where this increased complexity for atol is demonstrably better than the simpler atol::Real=0 because there is a loss of readability changing 0 to zero(promote_type(real(promote_leaf_eltypes(x)),real(promote_leaf_eltypes(y)))).

[thchr]

Okay, we can keep the original atol: I don't feel strongly either way.

The current asymmetry between the care taken with the types of rtol and atol is odd though: in one case, we try quite hard, in the other, not at all.
Incidentally, the explicit ::Real type-constraints on rtol and atol seem unfortunate: e.g., this - and the atol = 0 default makes it so that dimensionful numbers (e.g., Furlong or Unitful.jl values) cannot fall back to Base's isapprox methods.
In fact, Unitful.jl seems to currently copy-extend this entire method; I assume partly for these reasons (e.g., because a unitful number cannot compare to a non-unitful one; and also because an AbstractQuantity isn't a subtype of Real).

But those details are orthogonal to the main point of this PR, so let's just go with the original atol = 0.

[LilithHafner]

Improving the default atol to obviate that definition in Unitful.jl would be a nice contribution in a different PR.

This one looks good to me as it. Thanks!

[LilithHafner]

CI failures look unrelated