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 (JuliaLang/julia#23194, JuliaLang/julia#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 JuliaLang/julia#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 #55 (comment). 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 (JuliaDiff/DualNumbers.jl#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 andyferris/Traitor.jl#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 JuliaLang/julia#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? #55 (comment)

[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.

[timholy]

@tkf, what's the barrier to just using IntervalArithmetic for your use case?

[tkf]

@timholy It's rather aesthetic one. I just prefer minimalistic dependency if possible and also conceptually I'm using interval-as-a-set rather than interval-as-an-approximation-of-number. I'm using IntervalSets.jl with a few custom operators and I'll keep doing so if this PR is not merged. I just thought I'd upstream my code in case it's useful for others. (Also, rather selfishly, I thought I could get some code review "for free".)

[dpsanders]

I view IntervalArithmetic.jl as exactly providing "calculations with sets", not "approximations of numbers".

Note that there is a recent package [NumberIntervals.jl] which throws errors (I believe) for operations like a < b for intervals that are not 0-width intervals that behave like numbers. cc @gwater

For 0..0.3 the point is that the real number 0.3 == 3//10 is not in this interval. We have been using .. to mean "the interval that the user really wanted, interpreting floating-point numbers as if they were true real numbers", giving the following:

julia> using IntervalArithmetic

julia> x = interval(0, 0.3)
[0, 0.3]

julia> 3//10 ∈ x
false

julia> x = 0..0.3
[0, 0.300001]

julia> @format full
Display parameters:
- format: full
- decorations: false
- significant figures: 6

julia> x
Interval(0.0, 0.30000000000000004)

julia> 3//10 ∈ x
true

I have been doubting recently if this is a good idea, since who knows what "the user really wanted". But I certainly believe that there should be some way of doing this using a representation like 0.3 and not requiring the user to think that this is actually 3//10.

[tkf]

Please see #55 (comment). The reason why I think it's doing interval-as-an-approximation-of-number is that IntervalArithmetic.jl enforces "rounding outward":

julia> x = @interval(-5, 5)
[-5, 5]

julia> sign(x)
[-1, 1]

Of course, it depends on the definition of calculations with sets. I'm using the definition that works with arbitrary sets, not just intervals:

c in f.(x, y) if and only if there exist a in x and b in y such that c = f(a, b)

The answer of sign(x) above would be Set([-1, 0, 1]) based on this definition. (Not that I'm arguing this is useful.)

[dpsanders]

I agree that the answer is not "correct" for true set operations, but I wouldn't say that's approximating a number (which number is - 1..1 approximating?)

Rather, it's just returning the smallest (in the best case) interval that contains all possible results. In this way you can keep everything within the same universe of "calculations with intervals".

[tkf]

I used the phrase "approximation of number" as I thought that was the main application of IntervalArithmetic.jl. I agree that "calculations with interval" is the best short description.

[timholy]

I gave this a review since you asked. It's beautiful work, as usual!

Not certain we want to merge this, but no harm in looking.

[timholy]

throw?

[timholy]

Perhaps we want

coretype(::Interval{L,R}) where {L,R} = Interval{L,R}

Then this could just be

broadcasted(::typeof(/), i::AbstractInterval, x) = coretype(i)(endpoints(i) ./ x)...)

[tkf]

@timholy Wow, I didn't expect that coming. Thanks a lot!

[Kolaru]

I think we can have with the same .. definition as in here in IntervalArithmetic. The trick to "correctness without guessing" is probably to intercept floating point literals before they parsed to have control on their rounding.

This is actually what we do with the @interval macro, escaping float literals to string and then parsing them with the correct rounding.

Regarding the problem of Interval <: Real, numerous discussions in IntervalArithemetic showed many shortcomings of this approach. In particular it is not clear what constitute the core properties of Reals.

Regarding for example ordering (a property one generally expect from real Reals), we can ask if comparisons (<, >, etc.) of Real should always return a boolean. And if yes, does being Reals implies things like (a < b) == !(a >= b)? Having the first but not the last property will make any code of the form

if a < b
    # Do something
else
    # Do something else
end

likely to silently and unexpectedly fail.

[gwater]

And if yes, does being Reals implies things like (a < b) == !(a >= b)? Having the first but not the last property will make any code of the form

if a < b
    # Do something
else
    # Do something else
end

likely to silently and unexpectedly fail.

Ideally ambiguous cases should raise exceptions in comparisons, see also the behavior Boost implements, and what I adapted for NumberIntervals.jl. In any case - if this package is meant for sets, why not use set comparisons, like precedes()?

[tkf]

Using broadcasting for set arithmetics is still under discussion JuliaLang/julia#33282. But, if we can do this, a neat way to write this is all(a .< b). This code works with Reals, too.