Define interval arithmetics using broadcasting

[tkf]

close #32

[timholy]

So I foresee at least one headache. There's surely an easier demo, but here's one I happened to know:

julia> x = 0..3
0..3

julia> 30 in ((((((((x .+ x) .+ x) .+ x) .+ x) .+ x) .+ x) .+ x) .+ x) .+ x
true

julia> y = x ./ 10
0.0..0.3

julia> 3 in ((((((((y .+ y) .+ y) .+ y) .+ y) .+ y) .+ y) .+ y) .+ y) .+ y
false

julia> 0.3 in y
true

Lack of associativity with floating-point arithmetic means that some things that should be axiomatic are violated. This is part of why https://github.com/JuliaIntervals/IntervalArithmetic.jl goes to pains to use BigFloats, but of course that also makes that package painfully slow for many applications.

This package was originally constructed to be the "non-controversial" core operations involving intervals. Given the problem above, we would have to ask whether introducing this behavior changes that guarantee. Maybe that's OK, but we shouldn't do it lightly.

CC @dpsanders

[dlfivefifty]

BroadcastStyle could be used here: FastIntervalStyle and RigorousIntervalStyle.

[tkf]

@timholy Isn't it a problem only if something becomes floating-point when you didn't ask? That is to say, if you already have 0.0 .. 0.3 it's probably OK if the edges become inaccurate?

(0 .. 3) ./ 10 is a bit of an edge case but how about asking people to use (0 .. 3) .// 10 in this case? Following the analogy that "0 .. 3 is a container of Ints", it feels natural that the broadcasted operations follows the promotion rule of /(::Int, ::Int).

BroadcastStyle could be used here: FastIntervalStyle and RigorousIntervalStyle.

@dlfivefifty Alternatively, how about defining broadcasted(::typeof(big), ::AbstractInterval) so that big.(x) ./ 10 works?

[timholy]

@timholy Isn't it a problem only if something becomes floating-point when you didn't ask? That is to say, if you already have 0.0 .. 0.3 it's probably OK if the edges become inaccurate?

I'm a bit sensitive to this because I ended up spending about a month of my life immersed in the joys of endpoints for ranges (https://github.com/JuliaLang/julia/pull/23194, https://github.com/JuliaLang/julia/pull/18777 among others). I had the same perspective you do, but people didn't like it not being exact (bugs were reported when they were not). That said, given the TwicePrecision stuff now in Base this would presumably be rather easier to adopt here.

[dpsanders]

Here's a simpler example:

julia> 0.9 ∈ (y .+ y .+ y)
false

In fact in IntervalArithmetic.jl we have used Float64 for a long time (unless you explicitly request BigFloat), and a while ago we worked hard to also make that fast, but still with correct rounding (using ErrorfreeArithmetic.jl).

Unfortunately the sum of 10 intervals that @timholy used above still takes 85 ns, compared to 3.5ns for this PR (and 1.8 μs for BigFloat intervals -- 20x slower). There doesn't seem to be any reasonable way around the overhead, unfortunately, without messing with the processor's rounding mode, which has its own problems and overhead.

Given that there is already a correct implementation available elsewhere, maybe it makes sense to have this non-guaranteed implementation here? But then you are never sure if your result is really correct, which is part / most of the point of interval arithmetic.

I do like the broadcasting syntax for interval-as-container.

[dpsanders]

Note that

julia> 3//10 ∈ y
false

[tkf]

I'm a total noob when it comes to rigorous interval arithmetic. But is it conceivable to prefer rounding in the opposite direction than IntervalArithmetic.jl? That is to say, is there a situation where you want that the result is guaranteed to not contain incorrect result?

I'm asking this because, if you take the definition c in f.(x, y) if and only if there exist a in x and b in y such that c = f(a, b) very literary, "rounding outward" also seems to be mathematically incorrect. (I'm just imagining "rounding outward" is what IntervalArithmetic.jl is doing. Sorry if I'm wrong.)

But then you are never sure if your result is really correct, which is part / most of the point of interval arithmetic.

My use case is just shifting around two numbers in floating-point precision so I don't need rigorous arithmetic personally. But this does not appear in a tight loop so I don't mind to use more involved definition either.

[timholy]

But is it conceivable to prefer rounding in the opposite direction than IntervalArithmetic.jl? That is to say, is there a situation where you want that the result is guaranteed to not contain incorrect result?

I don't think there's a useful situation in which that's true. Consider x .- x when x = 1..2. The result -1..1 is clearly conservative; if your intention is actually y - y for y in x, the actual answer would be 0..0. -1..1 is the "most conservative" answer (complete lack of correlation between the elements selected from x), and in such cases quibbling over an eps at the edge is not worth giving up the guarantee that you know the answer lies within the interval no matter how you interpret the underlying meaning of that operation.

[timholy]

Let me also make it clear that I do think there is room for a fast-and-loose interval arithmetic package that takes floating point at face value; I just think that a package that promises to "get the topology right, and leave it at that" is not the place for such a thing.

[dpsanders]

Actually there is a (not very well documented and probably not completely implemented nor tested) "no rounding at all" mode in IntervalArithmetic.jl that satisfies the requirement for a "fast and loose mode".

[dlfivefifty]

Note there is really a need for two different types:

1. Interval <: Real, that is interval arithmetic. Here we want rigorous bounds, and support for standard arthimetic which is pessimistic: (0..1) + (2..3) == (2..4) and sign(-2..2) == -1..1. Broadcasting will also behave like interval arithmetic, but this is a consequence of behaving like scalars, e.g., exp.(1..5) gives the same thing as exp(1..5) as this is a generalisation of exp.(5) == exp(5).
2. Interval <: Domain, in which case an interval behaves like a set. Here we may not want standard arithmetic and definitely do not want exp(1..5) to be defined. Here we do not want pessimistic bounds: sign.(-2..2) should not be -1..1 as we want to interpret it as {sign(x) for x in -2..2}, so should actually be Set([-2,0,2]). An example usage would be using exp.(im.*(0..π)) to represent a half-circle.

For (2) I would say standard floating point results are found. If someone wants rigorous interval arithmetic they should be using (1) as in IntervalArithmetic.jl.

That said, there would be no problem in having this broadcast behaviour implemented outside this package, say a new package IntervalSetsBroadcasting.jl or even just putting it in DomainSets.jl

[tkf]

Here we do not want pessimistic bounds

My point on the emphasis on the "only if" part is in line with this argument.

definitely do not want exp(1..5) to be defined

I suppose you mean that treating Interval as a container is the right choice? I actually started doubting this argument; or, rather, wondering what the semantics of broadcasting should be in general. My mental model is that it "joins" indices/keys of the containers, including arrays (and do something more with singleton axes). This seems to be a popular semantics, looking at https://github.com/JuliaLang/julia/issues/25904. However, the semantics here is taking a product of two sets and then apply the function on it (i.e., set1 .+ set2 == Set([x + y for x in set1, y in set2])). This deviation can be confusing as @timholy brought up in https://github.com/JuliaMath/IntervalSets.jl/pull/55#issuecomment-531610861. I see two paths:

1. Treat Interval as a non-container and non-Real abstract object; i.e., define + etc.
2. Generalize broadcasting semantics somehow to allow taking product.

Note that the second path involves fixing the behavior in Base:

julia> Set([1,2]) .+ Set([3,4])
2-element Array{Int64,1}:
 6
 4

[dlfivefifty]

The behaviour of broadcasting with Set is clearly ill-defined as it’s picking an arbitrary ordering. So perhaps it makes sense to first do a PR to Base to try to introduce a SetStyle broadcasting style to establish the convention.

[timholy]

@dlfivefifty,

Interval <: Real, that is interval arithmetic.

I'd rather this weren't true in the long run, and we reserved Real for things that really are. Intervals don't satisfy several of the properties that other real numbers satisfy. I think Interval <: Real is pragmatically necessary right now, but once we have a good traits system I hope we can get rid of that.

@dpsanders, I learned some stuff about IntervalArithmetic! I'll have to get back to playing with it at some point, it sounds like some of the things I was concerned about are probably not an issue any more.

[dlfivefifty]

This "is it a Real" issue also comes up with dual numbers (https://github.com/JuliaDiff/DualNumbers.jl/issues/45). The reason Interval <: Real makes sense is that the usage of interval arithmetic often consists of applying functions with definitions of the form

f(x::Complex) = ...
f(x::Real) = ...

Another example is one gets complex rectangles for free if Interval <: Real: complex(0..1,2..3) would be supported and a lot of features will automatically work as implementations reduce to real number operations.

Traits may provide an alternative solution in Julia v2.0 or so, but at the moment that is speculative. In the meantime, I think it's best not to think of mathematical definitions but of "interfaces". Interval <: Real makes sense if intervals implements the (unwritten) real number interface.

[timholy]

I think it's best not to think of mathematical definitions but of "interfaces". Interval <: Real makes sense if intervals implements the (unwritten) real number interface.

I am not so sure you can productively separate these. If you write a bit of code thinking about "real Real" numbers, then the following will likely surprise you:

julia> x, y = 1..2, 1.5..1.8
([1, 2], [1.5, 1.80001])

julia> x < y, y < x, x == y
(false, false, false)

There are no MethodErrors here (the interface is defined), yet the consequence is likely to be a source of bugs that is much harder to diagnose.

As far as getting traits, I think someone just needs to spend a month on https://github.com/andyferris/Traitor.jl/issues/8. That doesn't help any definitions in Base, unfortunately, but I don't think there are any technical barriers for the package ecosystem to develop a nice trait system.

[dlfivefifty]

The problem is you often need to use interval arithmetic for other package codes. Forcing all packages to use traits to describe real numbers is not realistic, and constructions like Complex{<:Interval} will still be banned.

In any case, whether its a good idea or not, it's already being done:

https://github.com/JuliaIntervals/IntervalArithmetic.jl/blob/1516c5e1ca8e088aeafa9574072512aa2df241fb/src/intervals/intervals.jl#L14

I think perhaps the long term solution is to make interval an interface, not a type. So perhaps there could be IntervalBase.jl with a few basic functions that other packages implement, and then IntervalSets.jl will be the home of "intervals as sets" and IntervalArithmetic.jl will be the home of "intervals as numbers", with a shared interface:

module IntervalBase
"A tuple containing the left and right endpoints of the interval."
function endpoints end

"The left endpoint of the interval."
leftendpoint(d) = endpoints(d)[1]

"The right endpoint of the interval."
rightendpoint(d) = endpoints(d)[2]

"A tuple of `Bool`'s encoding whether the left/right endpoints are closed."
function closedendpoints end

"Is the interval closed at the left endpoint?"
isleftclosed(d) = closedendpoints(d)[1]

"Is the interval closed at the right endpoint?"
isrightclosed(d) = closedendpoints(d)[2]

# open_left and open_right are implemented in terms of closed_* above, so those
# are the only ones that should be implemented for specific intervals
"Is the interval open at the left endpoint?"
isleftopen(d) = !isleftclosed(d)

"Is the interval open at the right endpoint?"
isrightopen(d) = !isrightclosed(d)

# Only closed if closed at both endpoints, and similar for open
isclosed(d) = isleftclosed(d) && isrightclosed(d)
isopen(d) = isleftopen(d) && isrightopen(d)

function infimum end
function supremum end
end

[timholy]

In any case, whether its a good idea or not, it's already being done:

I know. It's true of ForwardDiff too. I'm just speaking about where we want to go, and less concerned about "what we have now."

So perhaps there could be IntervalBase.jl with a few basic functions that other packages implement, and then IntervalSets.jl will be the home of "intervals as sets"

In a sense I was originally envisioning that IntervalSets could be the base. IntervalArithmetic defines its own intersect, union, and similar, and I have a hard time envisioning any interval package that wouldn't need topology/set operations. That was part of what I had in mind with encouraging type-piracy. And also why I think one needs to be a bit careful about growing this package in directions that might interfere with some of its potential universality.

[timholy]

WRT to the original question of type piracy (@tkf), let me point out: having another package that "pirates" this one to add arithmetic operations would be a lot like ColorVectorSpace pirating the types defined in ColorTypes and adding arithmetic. I can't remember a single problem that's ever caused.

[dlfivefifty]

** In a sense I was originally envisioning that IntervalSets**

Yes I know. I think it didn't work out quite as well as anticipated as IntervalArithmetic.jl doesn't use it, so I'm proposing an even more "bare bones" package, just so there is a common interface (common type hierarchy appears to be impossible).

What's unclear is what to do about ... Perhaps there is no reason to make this shared so there would be both an IntervalSets: .. and IntervalArithmetic: .. and the downstream user can import the one they want.

The end result would be a package like ApproxFun.jl can mix-and-match IntervalSets.jl (for the function domains) and IntervalArithmetic.jl (for computations), using a common interface (endpoints, etc.) to query the relevant types, and do import IntervalSets: .. so that a..b always defaults to the function domain case.

[timholy]

Yes, a common interface sounds good in principle. I'm a bit unclear about what would go in it, but if it's practical...

[dlfivefifty]

Basically what I said above: endpoints, closedendpoints, etc. So very minimal. (There are problem others I forgot.)

[dpsanders]

In IntervalArithmetic.jl we (currently) use x..y to give an interval that is guaranteed to contain the real numbers x and y, e.g. 0..0.3 gives a result that is different from IntervalSets.jl. This is the main thing holding up using IntervalSets as a base package.

How problematic is it if we pirate .. to use this different (by up to one ulp) definition?

[tkf]

@dlfivefifty I opened https://github.com/JuliaLang/julia/issues/33282 asking if it makes sense to use broadcasting for set arithmetic. Can you comment there on why you think it is a good way to go, even though the connection to broadcasting with associative data structure is unclear (OK, at least for me)?

[tkf]

Until JuliaLang/julia#33282 is resolved, how about adding setadd, setmul etc. and some unicode operator shorcuts? They can be easily hooked into broadcasting mechanism once JuliaLang/julia#33282 is resolved.

The end result would be a package like ApproxFun.jl can mix-and-match IntervalSets.jl (for the function domains) and IntervalArithmetic.jl (for computations)

This would be great.

[tkf]

What's unclear is what to do about ...

If it were just .. (and ±?), how about just using the definition from IntervalArithmetic.jl in IntervalSets.jl? For fast (re)constructions, you can still call the constructor directly.

[dlfivefifty]

The issue is there are two different interval types and .. needs to return one or the other: one is <: Real and the other is <: Domain. We could in theory drop the latter, though then it’s propagating the lie that an interval is a real number.

[tkf]

Ah, yes, of course.

[tkf]

What do you think about adding custom functions/operators for now? https://github.com/JuliaMath/IntervalSets.jl/pull/55#issuecomment-531977627

[timholy]

@dpsanders,

we (currently) use x..y to give an interval that is guaranteed to contain the real numbers x and y, e.g. 0..0.3 gives a result that is different from IntervalSets.jl. This is the main thing holding up using IntervalSets as a base package.

Didn't know that, that's really interesting. Ours also contain the endpoints, of course, but I guess you're saying you "defensively" widen the interval a bit? I guess it comes down to whether you want "safe" or "literal" by default.

I bet the <:Real is a bigger obstacle, and one that might require better number traits to overcome.