Example 5.5 from Smiley and Chun (2001) solved inconsistently

[gwater]

The example in test/smiley_examples.jl has 41 roots, however, only 20 are found on some setups. On other setups all 41 are identified correctly. This inconsistency was first discovered and discussed in #70.

[gwater]

Some progress on triangulating the issue: On my machine the example works with v0.7.0 of StaticArrays.jl and breaks with v0.7.2.

[gwater]

The example is already broken with v0.7.1 of StaticArrays.jl

[gwater]

git bisection indicates that the problem is related to this commit in StaticArrays.jl

[dpsanders]

Great detective work, thanks. Would you like to report it at StaticArrays.jl?

[gwater]

Sure, it cannot hurt to have them take a look at it. The code in question is called here, correct?

[dpsanders]

Yes, that's the right place.

[gwater]

Given that the next version of StaticArrays will most likely contain a fix for this issue, I think v0.7.1 through v0.7.2 of StaticArrays.jl should be blacklisted in IntervalRootFinding.jl's REQUIRE to protect users from incorrect results.

[gwater]

The discussion in the StaticArrays issue raised questions about the definition of isinf(::Interval). Looking at its current definition I'm concerned, too:

isentire(x::Interval) = x.lo == -Inf && x.hi == Inf
isinf(x::Interval) = isentire(x)

This definition seems both geneerous and too specific at the same time. My understanding is that the question "Is this number infinite?" cannot be answered at all for an interval that contains both an infinity and finite numbers. It could be both finite or infinite, and we don't know yet.

Should I open an issue to discuss this further?

[dpsanders]

Yes, I agree. Please do open an issue, yes.

[dpsanders]

I think the solution for us is to use the work in #67. cc @eeshan9815

[gwater]

Yes, it makes sense to have a safe implementation of left divide within the package since that operation is so crucial to the Newton method.

[eeshan9815]

So should the methods in #67 be used instead of Base.\ by importing and modifying it?

[dpsanders]

Yes, that's probably a good solution. Let's try it.

[dpsanders]

Maybe make it do Gaussian elimination for now. But probably there should be a way to choose in the roots function which type of linear solution to use.

[eeshan9815]

Should it use the preconditioned version or the normal one?

[dpsanders]

There should be an option to either use it or not, and test both to see which is better!

[dpsanders]

I think this will be solved by using interval Gauss-Seidel instead of Gaussian elimination to solve the linear system. @Kolaru is it possible to specify this option?

[Kolaru]

@dpsanders Currently it's not possible, the contractors using the \ operator (and thus whatever it is bind to).

I've added it to the list in #116 of what should be possible to pass to roots.

[gwater]

Silent failure should be solved by JuliaIntervals/IntervalArithmetic.jl#571

[gwater]

With the earlier change to the contractors (using a custom \ operation) and the restriction on boolean operations in IntervalArithmetic.jl this issue can be considered resolved. Both precautions together should be sufficient.