overload isapprox for intervals

[lucaferranti]

isapprox wasn't overloaded for intervals and was using the fallback isapprox(a::Number, b::Number). This PR defines the method isapprox(a::Interval, b::Interval)

before

julia> a = 0.3..0.7
[0.299999, 0.700001]

julia> b = a + 1e-10
[0.3, 0.700001]

julia> a ≈ b
false

after

julia> a = 0.3..0.7
[0.299999, 0.700001]

julia> b = a + 1e-10
[0.3, 0.700001]

julia> a ≈ b
true

Related issues:

fixes #479

[codecov-commenter]

Codecov Report

Merging #480 (1ebe004) into master (c7602f4) will decrease coverage by 0.00%. The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #480      +/-   ##
==========================================
- Coverage   91.56%   91.55%   -0.01%     
==========================================
  Files          25       25              
  Lines        1755     1753       -2     
==========================================
- Hits         1607     1605       -2     
  Misses        148      148              

Impacted Files	Coverage Δ
src/IntervalArithmetic.jl	100.00% <ø> (ø)
src/intervals/arithmetic.jl	97.21% <100.00%> (+0.01%)	:arrow_up:
src/multidim/intervalbox.jl	88.88% <0.00%> (-0.59%)	:arrow_down:

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update c7602f4...1ebe004. Read the comment docs.

[Kolaru]

This is one of the case where we need to be very careful. There are two possible strategy (already discussed several times in e.g. #181):

(a) We do the most straightforward thing, so for this case two intervals are approximately equals if their respective bounds are. This is what is implemented here. However it may break code written for numbers in unexpected ways. (b) We extend all functions defined on numbers in the interval sense so that to ensure that if you replace a number by an interval you get a guaranteed correct result. In this case a ≈ b should probably be missing.

This is also related to https://github.com/gwater/NumberIntervals.jl as an actual implementation of (b) and #271 that tried to build the ground work to solve all these issues. Unfortunately I don't know when I will have the time/energy to go back to the latter.

[lucaferranti]

Hups, apologies for rushing into the PR without the proper background check, I was aware of the issues/PRs you linked but when I opened this I did not connect the dots.

I can close this PR for the time being, but I do have a couple of comments:

it may break code written for numbers in unexpected ways.

This is happening with the current version, which is falling back to Base.isapprox, which uses norm which for objects of type Number uses the absolute value, which for intervals is an interval, which makes no sense in this contest.

Regardless of what the correct behavior of isapprox should be, the current one is not,

julia> a = 1..(1+1e-16)
[1, 1]

julia> b = a + 1e-308
[1, 1.00001]

julia> a ≈ b
false

julia> a.lo ≈ b.hi
true

I chose this example because in this case any number in a is approximately equal to any number in b (with the default tolerance settings), so this should be true in both philosophies. Always returning missing for == and isapprox would cut the head off the snake, but I am a little skeptical it would be practically useful in bigger applications. A proper implementation of (b) would require modal logic I think, something like [1, 2] ≈ [1, 2] is possibly true and a ≈ b from the example above is definitely true.

Moreover, I would say that generic code written for Number would break not because some functions are overloaded for Interval as here using philosophy (a), but because some functions aren't, but if you have examples in mind where overloading == and isapprox breaks, please do prove me wrong :)

[lucaferranti]

Also, reading at the old related issues, my main take-home understanding (but I wouldn't be surprised if I misunderstood something) is that:

while the package is mainly following philosphy (a), having Interval <: Real is a contradiction, because intervals thought as sets do not fulfill the axioms of real numbers`.

I agree with this! Technically speaking, floating-point numbers do not obey the axioms of real numbers either, since e.g. NaN has no additive inverse. However, it is true that not even "ideal interval" with infinite precision etc. obey the axioms of real numbers e.g. distributivity does not hold.

Hence I agree that having Interval <: Real leads to a semantical contradiction and I understand why considering philosophy (b).

However, what about having Interval <: Number I am not an expert in type/set theory or foundations of mathematics in general, but I would say that numbers (whatever a number is) do not have any particularly strict axiom that would be broken by intervals and operations as defined here following (a)?

I am sure you had already thought about this and I apologize if it feels like we are having just another instance of an old conversation. I will re-read the discussions in the related issues and PRs to make sure I don't waste anyone's time by reinventing the wheel. For what is worth, I do think that repeating an old conversation after some time is not necessarily a waste of time, as it can bring up new insights 😄

[Kolaru]

First no need to apologize, I pointed out old issues because I thought you may find them interesting, not because I expected you to know about them :) (this goes as far back as #2 so it is quite a rabbit hole)

Always returning missing for == and isapprox would cut the head off the snake, but I am a little skeptical it would be practically useful in bigger applications.

I think that underlines the core of the problem. You are right, always returning missing is not useful. It will simply break any code that expects == or isapprox to return a boolean. But this is safe. If the code never runs, it will never produce wrong result. The tricky part is that we do not have any other way to guarantee that the code that was given Interval is producing guaranteed result when boolean comparison are involved.

if you have examples in mind where overloading == and isapprox breaks, please do prove me wrong :)

An example would be a possible Taylor expansion of a function f(a -b)

if a ≈ b
    # Do a Taylor expansion of the function around x
    return f(0) + (a - b) f'(0)
else
    # Go the expensive route
    return f(a - b)
end

If you input the interval in you example a = b = (3..7) the first branch will be used, but a - b == (-4..14) which can not be considered to be close to zero and is certainly not proper for a first order Taylor expansion.

It relies of one of the properties of numbers that intervals are missing (a == b => a - b == 0).

That kind of branching is not unusual for the computation of special functions like the zeta function (however it generally uses inequalities I think).

what about having Interval <: Number

A lot of packages require <:Real to signal they want a number that is not a complex number. Removing Interval <: Real would prevent interoperability. Currently the consensus (IIRC) in #271 was to implement strict rules for boolean functions, at least for the default flavor, while keeping Interval <: Real.

Another course of action (briefly discussed on slack) would be to encourage removal of this type restriction everywhere and mostly rely on duck typing. But that would be a more global discussion. This is however tangential to the problem of boolean comparison of intervals.

[lucaferranti]

thank you for the example, it was an interesting one, also apologies for returning to this so seldomly. I promise that next time I'll answer faster 😄

The current implementation of == at the moment is

function ==(a::Interval, b::Interval)
    isempty(a) && isempty(b) && return true
    a.lo == b.lo && a.hi == b.hi
end

which is in accordance to philosophy a) of your original message. This may change in the future and if in the future we will have the flavour based system with different mechanisms, but since at the moment the "set flavor" is being used for equality, I don't understand why not use the same for approximate inequality, that is here in the PR I am mainly after coherence in the package.

The current fall back implementation is failing because of

norm(x-y) <= max(atol, rtol*max(norm(x), norm(y))))

and norm(a::Interval) is just returning abs(a), which is an interval and not a positive real.

[lucaferranti]

Also, regarding packages interoperability, preventing it is not always a bad thing. For example, plugging interval matrices into traditional linear algebra algorithms is usually a bad idea, let's take an example

julia> using IntervalArithmetic, LinearAlgebra

julia> A = [1..2 1..2;2..2 3..3]
2×2 Matrix{Interval{Float64}}:
 [1, 2]  [1, 2]
 [2, 2]  [3, 3]

julia> lu(A)
LU{Interval{Float64}, Matrix{Interval{Float64}}}
L factor:
2×2 Matrix{Interval{Float64}}:
 [1, 1]  [0, 0]
 [1, 2]  [1, 1]
U factor:
2×2 Matrix{Interval{Float64}}:
 [1, 2]   [1, 2]
 [0, 0]  [-1, 2]

everything seems nice we love having packages just magically work. But let us take now a nastier example

julia> A = [1e-16..2 1..2;2..2 3..3]
2×2 Matrix{Interval{Float64}}:
     [1e-16, 2]  [1, 2]
 [2, 2]          [3, 3]

julia> lu(A)
LU{Interval{Float64}, Matrix{Interval{Float64}}}
L factor:
2×2 Matrix{Interval{Float64}}:
 [1, 1]       [0, 0]
  [1, 2e+16]  [1, 1]
U factor:
2×2 Matrix{Interval{Float64}}:
     [1e-16, 2]  [1, 2]
 [0, 0]             [-4e+16, 2]

[-4e16, 2]? that's not a nice interval! So what happened? If we @edit lu(A) we see that the generic implementation uses the maximum absolute value as pivot, because you want to choose as pivot a number as far as possible from zero, and the as far as possible is expressed by maximum absolute value for real numbers, but for intervals you cannot use maximum absolute value, but maximum mignitude. In the example below, you initialize the first column pivot to [1e-16, 2] and then you check abs([2, 2]) >= abs([1e-16, 2]), which evaluates to false and hence you use [1e-16, 2] as pivot, divide by a number close to zero and everything goes bananas. (btw, proper gaussian elimination with maximum mignitude pivoting is implemented in IntervalLinearAlgebra.jl)

This and the fact that <: Real was mainly introduced for interoperability with ForwardDiff and it's not required anymore makes me think that <:Real should be removed (or replaced to <: Number), but yeah I think that's tangential to this PR

[OlivierHnt]

@lucaferranti what is the status of this PR? Perhaps rebase and define some isapprox_interval function?

[OlivierHnt]

I am closing this PR since it seems abandoned. Feel free to re-open it if you want to revive it at some point, or just open a new one.