Root finding with interval parameters

[jebej]

I couldn't find an example for this in the docs, but are Interval values (not unknowns) allowed in the root function?

I tried the following (note the interval g1 parameter), but it doesn't appear to finish. Any advice on if this is doable?

Note that the code works great if I plug in a number for g2, but I would like the result interval to reflect the uncertainty in g2.

using IntervalRootFinding, IntervalArithmetic, StaticArrays

function f(z)
    # units in MHz
    f1 = 5303.5  # frequency for resonator 1
    χ1 = 0.34144 # shift for resonator 1
    g1 = 38..42      # estimated coupling for resonator 1

    f2 = 4522.7  # frequency for resonator 2
    χ2 = 0.24841 # shift for resonator 2

    fq, g2 = z   # unknown qubit freq and coupling for resonator 2
    return SVector(
        (fq-f1)*χ1 + g1^2,
        (fq-f2)*χ2 + g2^2
    )
end

rts = roots(f, IntervalBox(0..8000, 0..100))

[dpsanders]

In principle yes they are allowed, but the problem is that the zeros are then no longer isolated, so the standard interval Newton method (which roots uses by default) will not work.

Using bisection instead should work:

rts = roots(f, IntervalBox(0..8000, 0..100), Bisection)

I believe that this could be solved by using the parametric interval Newton method, which should be able to prove uniqueness for each parameter value in your input interval, but that is not yet implemented. Help here is welcome!

[Kolaru]

It seems to run infinitely because the default tolerance is absolute and small (1e-7). If you put in a tolerance more adapted to your range it terminates

julia> roots(f, IntervalBox(0..1000, 0..100), Newton, 10.0)
111-element Array{Root{IntervalBox{2,Float64}},1}:
 Root([792.702, 800.575] × [28.9393, 33.6693], :unknown)
 Root([299.76, 307.633] × [28.9393, 33.6693], :unknown)
 Root([338.634, 346.385] × [28.9393, 33.6693], :unknown)
 Root([330.762, 338.635] × [28.9393, 33.6693], :unknown)
 Root([511.235, 518.866] × [28.9393, 33.6693], :unknown)
 Root([879.797, 887.67] × [28.9393, 33.6693], :unknown)
 Root([855.929, 863.926] × [28.9393, 33.6693], :unknown)
 Root([959.388, 967.511] × [28.9393, 33.6693], :unknown)
 Root([927.401, 935.524] × [28.9393, 33.6693], :unknown)
 Root([416.749, 424.622] × [28.9393, 33.6693], :unknown)
 Root([385.504, 393.255] × [28.9393, 33.6693], :unknown)
 Root([651.111, 658.862] × [28.9393, 33.6693], :unknown)
 ⋮
 Root([534.245, 541.996] × [28.9393, 33.6693], :unknown)
 Root([729.96, 737.957] × [28.9393, 33.6693], :unknown)
 Root([307.632, 315.263] × [28.9393, 33.6693], :unknown)
 Root([323.012, 330.763] × [28.9393, 33.6693], :unknown)
 Root([951.392, 959.389] × [28.9393, 33.6693], :unknown)
 Root([214.618, 222.369] × [28.9393, 33.6693], :unknown)
 Root([635.489, 643.24] × [28.9393, 33.6693], :unknown)
 Root([175.744, 183.617] × [28.9393, 33.6693], :unknown)
 Root([238.112, 246.11] × [28.9393, 33.6693], :unknown)
 Root([253.62, 261.251] × [28.9393, 33.6693], :unknown)
 Root([604.241, 612.114] × [28.9393, 33.6693], :unknown)
 Root([919.405, 927.402] × [28.9393, 33.6693], :unknown)

[jebej]

Thanks! Yes I have been experimenting with the tolerance.

In this case, it seems strange to return all of these roots. I know that there will be a big range with one root, so I just defined:

Base.union(a::Root, b::Root) = union(a.interval,b.interval)
Base.union(a::IntervalBox, b::Root) = union(a,b.interval)

to get:

julia> reduce(union, roots(f, IntervalBox(0..8000, 0..100), Newton, 1.0))
[137.145, 1074.36] × [29.2653, 33.008]

julia> reduce(union, roots(f, IntervalBox(0..8000, 0..100), Bisection, 1.0))
[136.801, 1074.69] × [29.201, 33.0763]```

edit: I suppose mapreduce is the better way of doing this:

mapreduce(r->r.interval, union, roots(f, IntervalBox(0..8000, 0..100), Newton, 1.0))

[jebej]

I am fairly confident that all the roots returned here overlap, so should the answer there be just a single root?

E.g. with a bigger tol so that the list is more manageable:

julia> roots(f, IntervalBox(0..1000, 0..100), Newton, 100.0)
14-element Array{Root{IntervalBox{2,Float64}},1}:
 Root([496.093, 557.618] × [0, 42.7511], :unknown)
 Root([808.57, 872.048] × [0, 42.7511], :unknown)
 Root([935.523, 1000] × [0, 42.7511], :unknown)
 Root([307.632, 370.125] × [0, 42.7511], :unknown)
 Root([746.078, 808.571] × [0, 42.7511], :unknown)
 Root([183.616, 246.11] × [0, 42.7511], :unknown)
 Root([557.617, 620.11] × [0, 42.7511], :unknown)
 Root([872.047, 935.524] × [0, 42.7511], :unknown)
 Root([370.124, 432.618] × [0, 42.7511], :unknown)
 Root([246.109, 307.633] × [0, 42.7511], :unknown)
 Root([137.145, 183.617] × [0, 42.7511], :unknown)
 Root([432.617, 496.094] × [0, 42.7511], :unknown)
 Root([620.109, 682.602] × [0, 42.7511], :unknown)
 Root([682.601, 746.079] × [0, 42.7511], :unknown)

Why are these intervals even being subdivided in the first place?

[Kolaru]

Interval bodes are being subdivised until one of the following is true:

• The radius of the box is smaller than the tolerance.
• The box is proved to not contain any solution (the box is discarded).
• The box is proved to contain a unique solution (the box is refined and saved).

The core idea is that proving existence or absence of solution is easier with smaller box.

Reassambling boxes is interesting, but it is not clear how it should be done in general. Taking the union results in a bigger box than necessary in many case.

Finally another issue here is that we only allow the same tolerance for each dimension, which in your case doesn't quite make sense.

[jebej]

Thanks for the detailed answer.

Interval bodes are being subdivised until one of the following is true:

• The radius of the box is smaller than the tolerance.
• The box is proved to not contain any solution (the box is discarded).
• The box is proved to contain a unique solution (the box is refined and saved).

The core idea is that proving existence or absence of solution is easier with smaller box.

This makes sense!

Reassambling boxes is interesting, but it is not clear how it should be done in general. Taking the union results in a bigger box than necessary in many case.

I suppose that as long as the intervals are overlapping, as they are for the above, it is natural to merge them, no? There could be an additional pass at the end of the root finding algorithm to do that. (I assume that when splitting a box, the resulting subboxes are created overlapping on purpose?).

I could look into doing that, is there already functionality to check if boxes are overlapping?

Finally another issue here is that we only allow the same tolerance for each dimension, which in your case doesn't quite make sense.

Yes, that could be solved by allowing to pass a tuple for the tolerance in each dimension.

[Kolaru]

I suppose that as long as the intervals are overlapping, as they are for the above, it is natural to merge them, no?

It totally is, I am just cautious before saying it will be straightforward :). Consider the following image where gray boxes have status :unkown and red boxes are known to be empty.

How should the gray boxes be merged? Taking the union would give big box covering everything, thus losing resolution, which is not acceptable (by default at least). Still you can merge 1-3 and 4-6, but you could also merge 3-4-6.

The most reasonable think to do may be to walk the internal tree in reverse and merge leaves with common ancestors. But really before someone try to seriously tackle that it is impossible to tell.

If you want to be the one tackling it you are of course welcome :)

[dpsanders]

The boxes produced by the algorithm don't really overlap, they abut: they share only edges (no interior points). They have to share edges in order to cover all possible inputs (since we deal only with closed intervals, not open intervals).

I wrote some code to unify boxes using connected components from LightGraphs.jl, which I can dig up, but as @Kolaru points out, the result is a single box that is the interval hull of the subboxes. Indeed it would be more sensible to use the tree to merge neighbouring boxes.

[dpsanders]

Ah in fact a version of unify is here: https://github.com/JuliaIntervals/IntervalOptimisation.jl/pull/5/files