Correct multidim Krawczyk

[dpsanders]

This "corrects" multi-dimensional Krawczyk to use the inverse of a real matrix instead of an interval matrix. This should not affect correctness but should significantly improve performance.

cc @Kolaru

[Kolaru]

This introduces a regression with respect to #93, since the center of the interval may be outside the definition domain of the Jacobian. The easiest way to deal with it is probably to use a try...catch statement and if a DomainError is caught return entireinterval. I am not a big fan of try...catch but the alternative is to compute the domain of the function beforehand and that looks hard to implement.

Also is there a reason to not adapt the 1D operator the same way ?

[dpsanders]

I thought that the 1D version already was like this, but I was wrong, so that should be done too.

OK I didn't think about the subtlety -- let me think about it. We definitely shouldn't use try.

[dpsanders]

What do you mean by the "definition domain"?

[Kolaru]

By "definition domain" I mean the domain of a function (it's a bad translation).

the definition domain of the Jacobian.

Precisely, here I mean the domain of the function jacobian(x) that gives Jacobi matrix of some function f at x. This domain is not the same as the domain of f in general.

[dpsanders]

You mean when f is not differentiable, e.g. for abs?

[Kolaru]

Yes for example. The example I had in mind was f(x) = x*log(x), which derivative (f'(x) = lox(x) + 1) has 0..Inf as domain.

[dpsanders]

And yet

julia> f(x) = x * log(x)
f (generic function with 1 method)

julia> f(-1..1)
[-∞, ∞]

julia> ForwardDiff.derivative(f, -1..1)
[-∞, ∞]

[dpsanders]

(I like the function lox(x) by the way ;) )

[dpsanders]

Hmm, this is nasty.

We in fact already have the machinery available for calculating the safe domain, which is used in the DecoratedIntervals stuff:

julia> f(DecoratedInterval.(X))
2-element SArray{Tuple{2},DecoratedInterval{Float64},1,2}:
     [-∞, ∞]_trv
 [-300, 300]_com

julia>     X = IntervalBox(-100..100, 2)
[-100, 100] × [-100, 100]

julia> f(xx) = ( (x, y) = xx; SVector(log(y/x) + 3x, y - 2x) )
f (generic function with 1 method)

julia> f(DecoratedInterval.(X))
2-element SArray{Tuple{2},DecoratedInterval{Float64},1,2}:
     [-∞, ∞]_trv
 [-300, 300]_com

Here, the trv decoration shows that it realises that there is "something wrong"; this came from calculating the domain. Maybe we can use this in a cleverer way.

[dpsanders]

I see why you suggested try...catch..., but from what I understand that is terribly slow.

[dpsanders]

Ah, how about this: Just check if f(IntervalBox(m)) is empty!

[Kolaru]

It can only work if we can be sure that the domain of the Jacobian is contained in the domain of the derivative. It is not the case for abs(x) for example. However this will only happen for non differentiable functions, and supporting piecewise differentiable functions seem to be enough for practical use (I doubt anybody will use numerics with crazy non smooth functions).

Dealing with these corner cases will need a bit of work though, your solution plus vibrant warning in the documentation seems sufficient for the time being.

[dpsanders]

@Kolaru How do you specify where to take the midpoint now?

[Kolaru]

Here: https://github.com/JuliaIntervals/IntervalRootFinding.jl/blob/dc41f395cb4939e0a7d58b855f6731f4508a8d27/src/contractors.jl#L123

It would make sense to make it go up to the top (i.e. roots function or the planned RootsProblem).

[dpsanders]

Thanks. The problem in this PR seems to be that it is bisecting exactly at 0 since the intervals are infinite.

On Thu, Mar 7, 2019 at 7:37 PM Benoît Richard notifications@github.com wrote:

Here: https://github.com/JuliaIntervals/IntervalRootFinding.jl/blob/dc41f395cb4939e0a7d58b855f6731f4508a8d27/src/contractors.jl#L123

It would make sense to make it go up to the top (i.e. roots function or the planned RootsProblem).

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaIntervals/IntervalRootFinding.jl/pull/114#issuecomment-470770404, or mute the thread https://github.com/notifications/unsubscribe-auth/AALtTlN8jx5h6xFXKnQTv9hpPt51nphlks5vUb7jgaJpZM4baz7p .

-- Dr. David P. Sanders

Profesor Titular "B" / Associate Professor Departamento de Física, Facultad de Ciencias Universidad Nacional Autónoma de México (UNAM)

dpsanders@g dpsanders@ciencias.unam.mxmail.com / Twitter: @DavidPSanders https://twitter.com/DavidPSanders http://sistemas.fciencias.unam.mx/~dsanders / GitHub: dpsanders https://github.com/dpsanders

Cubículo / office: #414, 4o. piso del Depto. de Física

Tel.: (+52 55) 5622 4965

[dpsanders]

Krawczyk is failing with infinite intervals on x->x^2. I think we should make bisect with a second parameter (bisection position) not bisect at 0 for infinite intervals.

[Kolaru]

Well this is currently the case, bisect can take two arguments and do not always bisect infinite intervals at 0:

julia> bisect(-Inf..Inf, 0.3)
([-∞, -7.19077e+307], [-7.19078e+307, ∞])

I may have integrated it in IntervalRootsFinding.jl along the way in #97, because with master I got:

julia> roots(x -> x^2, -Inf..Inf, Krawczyk)
1-element Array{Root{Interval{Float64}},1}:
 Root([-2.92648e-17, 8.14937e-16], :unknown)

The fact that the status is unknown is expected as the derivative at the solution is zero. If we remove that everything seems fine.

julia> roots(x -> x*(x - 1), -Inf..Inf, Krawczyk)
2-element Array{Root{Interval{Float64}},1}:
 Root([0.999999, 1.00001], :unique)
 Root([-1.72771e-18, 1.75491e-18], :unique)

So maybe rebasing into master will fix several of the issues.

[dpsanders]

I still get

julia> roots(x->x^2 - 2, -Inf..Inf, Krawczyk)
0-element Array{Root{Interval{Float64}},1}

Krawczyk uses mid, rather than bisect, which seems to be the problem.

[Kolaru]

In fact it has nothing to do with the bisection, that never land on exactly 0:

julia> f(x) = x^2
f (generic function with 1 method)

julia> contractor = Krawczyk(f, x -> 2*x)
Krawczyk{typeof(f),getfield(Main, Symbol("##3#4"))}(f, getfield(Main, Symbol("##3#4"))())

julia> search = BreadthFirstSearch(-Inf..Inf, contractor, 1e-10)
BreadthFirstSearch{Interval{Float64},Krawczyk{typeof(f),getfield(Main, Symbol("##3#4"))},Float64}(Root([-∞, ∞], :unknown), Krawczyk{typeof(f),getfield(Main, Symbol("##3#4"))}(f, getfield(Main, Symbol("##3#4"))()), 1.0e-10)

julia> global endtree

julia> for (k, tree) in enumerate(search)
           println("Iteration $k")
           println(tree)
           global endtree = tree
       end
Iteration 1
Working tree with 3 elements of type Root{Interval{Float64}}
Indices: [1, 2, 3]
Structure:
  [1] Node with children [2, 3]
    [2] Leaf (:working) with data Root([-∞, -1.40444e+306], :unknown)
    [3] Leaf (:working) with data Root([-1.40445e+306, ∞], :unknown)

Iteration 2
Working tree with 2 elements of type Root{Interval{Float64}}
Indices: [1, 3]
Structure:
  [1] Node with children [3]
    [3] Leaf (:working) with data Root([-1.40445e+306, ∞], :unknown)

Iteration 3
Working tree with 4 elements of type Root{Interval{Float64}}
Indices: [1, 3, 4, 5]
Structure:
  [1] Node with children [3]
    [3] Node with children [4, 5]
      [4] Leaf (:working) with data Root([-1.40445e+306, 8.84748e+307], :unknown)
      [5] Leaf (:working) with data Root([8.84747e+307, ∞], :unknown)

Iteration 4
Working tree with 3 elements of type Root{Interval{Float64}}
Indices: [1, 3, 5]
Structure:
  [1] Node with children [3]
    [3] Node with children [5]
      [5] Leaf (:working) with data Root([8.84747e+307, ∞], :unknown)

Iteration 5
Working tree with 0 elements of type Root{Interval{Float64}}

The problem is that the status of the interval [-1.40445e+306, 8.84748e+307] (iteration 3 on node 4) is incorrectly determined, probably because the contraction is broken:

julia> import IntervalRootFinding: 𝒦, where_bisect

julia> X = -1.40445e+306..8.84748e+307
[-1.40446e+306, 8.84749e+307]

julia> C = 𝒦(f, f′, X, where_bisect)
[-∞, -1.79769e+308]

julia> isempty(C ∩ X)
true

I think that the reason is that with that interval f(m) overflows to Inf. Adding the following check to the contractor seems to solve the issue:

val = f(m)
if isinf(val)
    return entireinterval(T)
end

[dpsanders]

Thanks for finding the bug -- and for a nice use of the tree introspection! I'll try adding that then.

[dpsanders]

@kolaru Does this version take account of your comments?

[Kolaru]

That seems good to me yes. Could you still add the problematic example we add to the tests?

This example was

roots(x->x^2 - 2, -Inf..Inf, Krawczyk)

that should find 2 unique roots.

[dpsanders]

Good idea, thanks. That does work; I'll add it as a test.