Incorrect unique zero for 0/x + x and similar.

[kemchurch]

I am computing zeroes of some pretty weird functions, and occasionally have unexpected behavior. A minimal working example (it's embedded in a lot of my functions; they depend on parameters and sometimes this kind of thing happens) is provided by the following.

julia> h(x) = x + zero(typeof(x))/x;

julia> r = roots(h,-1..1)
1-element Vector{Root{Interval{Float64}}}:
 Root([0, 0], :unique)

julia> h(r[1].interval)
∅

I think this is more of a feature of IntervalArithmetic.jl and/or ForwardDiff.jl internals than it is with IntervalRootFinding.jl, since the main reason this happens can be traced back to two things:

1. ForwardDiff.jl handles the derivative of zero(typeof(x))/x in such a way that zero(typeof(x))*ForwardDiff.derivative(y->1/y,x) is computed. This is just a hunch; I haven't looked at the internals.
2. IntervalArithmetic.jl handles certain products involving intervals with Inf and 0 in just the right way to cause trouble. Namely, zero(typeof(x))*interval(-Inf,Inf) will return [0,0], but interval(0)/interval(0) returns the empty interval. This is why interval Newton converges. Despite this, our function does not have a zero, and evaluates to the empty interval there as it should.

I suspect the same thing will happen if h is replaced with h(x) = zero(typeof(x))/x + f(x) for any f that has a unique zero in [-1,1] at 0 such that roots() converges.

I wasn't sure where to raise this issue. One could argue that this happens because of an inconsistency with IntervalArithmetic.jl's handling of 0*interval(-Inf,Inf) and similar. However, I noticed this problem while using IntervalRootFinding.jl, so I've raised the issue here for now.

[lbenet]

Hi @kemchurch! Thanks for opening this discussion here.

In general, I agree with your analysis. Yet, there is still one point you did not include stressed above, which i think is important: a naive evaluation of h(x) at 0.0 yields NaN. Clearly, the function $h(x)$ is discontinuous at $x=0$, and the discussion boils down to how is $h(0)$ defined; I am unsure that the interval Newton method theorem applies for this function. (Warwick Tucker, in his book "Validated Numerics: a short introduction to rigorous computations", proves the interval Newton method theorem (theorem 5.4) assuming that "$N([x])$ (Newton's operator) is well defined in the interval $[x]$".)

[kemchurch]

Ah, you make a good point. I suppose I was thinking of IntervalRootFinding.jl as something that it isn't: a black box that computes all zeroes of a given function on a box and encloses "bad" points (with status :unknown). It does indeed work like this in the case of some essential singularities (like for 1/x on a domain containing zero), but I shouldn't expect the same behaviour to hold generally.

[lbenet]

The idea of tagging roots as :unknown is because for the tolerance chosen they are not :unique (the Newton operator does not satisfy certain uniqueness property, but the existence of a root can not be ruled out), and may not be roots at all.

[Kolaru]

This main problem you are hitting here is that in the current flavor of IntervalArithmetic, $\infty$ is never part of an interval, since interval are subset of $\mathbb{R}$ and infinity is not a real number. Consequently [-Inf, Inf] is the open interval $(-\infty, \infty)$. So it is in fact consistent that 0 * [-Inf, Inf] == [0, 0], even so it is a bit weird since we still print it as a closed interval, and all other are closed intervals.

Allowing for different flavors of interval, and ultimately one with a more intuitive use of Inf (a candidate is the Cset type of interval), is already on the long list of things that would be cool to implement in IntervalArithmetic, but we are lacking the manpower for it.

I suppose I was thinking of IntervalRootFinding.jl as something that it isn't: a black box that computes all zeroes of a given function on a box and encloses "bad" points (with status :unknown).

In principle this could be true. We would need for it, as far as I can see:

• Using decorated interval to track things going out of the domain of definition of the function (this may already solve this particular issue by itself)
• The Cset type of interval mention above
• Proper handling of boolean comparison of intervals

So I think this issue should be kept open, as it is fair game to except it to work (or at least gracefully crash instead of giving wrong result) with the package.

[dpsanders]

As both @lbenet and @Kolaru have correctly pointed out, the Newton algorithm is designed to work with differentiable functions

Decorated intervals were designed to check whether a function over a given interval satisfies the requirements. (I thought we had some docs on this but I can't find them right now. cc @lucaferranti)

Here's your example use case. (Note that you can use zero(x) instead of zero(typeof(x)) and due to promotion you can actually just write 0/x!)

julia> using IntervalArithmetic

julia> @format true;   # display decorations

julia> f(x) = x + (0 / x);

julia> x = DecoratedInterval(1..2);

julia> f(x)
[1, 2]_com

julia> y = DecoratedInterval(-1..2);

julia> f(y)
[-1, 2]_trv

The decoration com ("common") shows that everything went fine, so the function f is continuous and bounded on 1..2. The decoration trv ("trivial") shows that "something went wrong" [or rather, technically, something may have gone wrong] and hence nothing can be trusted for the function f over -1..2. In this case it is clear that what went wrong is that 1/x is not defined for x==0, but for more complicated functions it may not be clear what went wrong.

[dpsanders]

The docs on decorations are here.

They don't seem to have been captured in juliaintervals.github.io. cc @lucaferranti