[Ann] IntervalRootFinding.jl for finding all roots of a multivariate function

Community
, ,
1

We are excited to announce v0.2 of the IntervalRootFinding.jl package, which has been completely rewritten from scratch.

This package exports a function roots that is guaranteed to find all roots of a given function \mathbb{R}^n \to \mathbb{R}^n in a given box X \subseteq \mathbb{R}^n (or tell you that it is unable to do so).
Functions are written with standard Julia syntax; no macros or other magic are required!

It is able to do this thanks to the power of interval arithmetic methods, based on the IntervalArithmetic.jl package. A generic “branch and prune” method is used to bisect space to look for roots.

As a flavour of what it can do, let’s find all the stationary points of the function f(x, y) = \sin(x) * \sin(y) in a box, by finding zeros of its gradient, i.e. vectors \mathbf{x} such that (\nabla f)(\mathbf{x}) = \mathbf{0}:


julia> using IntervalArithmetic, IntervalRootFinding

# define the function
julia> f(xx) = ( (x, y) = xx; sin(x) * sin(y) )

# obtain its gradient using the gradient operator ∇ exported by the package
julia> ∇f = ∇(f)

# calculate the roots of the gradient using the interval Newton method in the given box:
julia> rts = roots(∇f, (-5..6) × (-5..6), Newton, 1e-5)
25-element Array{IntervalRootFinding.Root{IntervalArithmetic.IntervalBox{2,Float64}},1}:
 Root([4.71238, 4.71239] × [4.71238, 4.71239], :unique)
 Root([4.71238, 4.71239] × [1.57079, 1.5708], :unique)
 ⋮
 [output snipped]

This returns a vector of pairs (interval, status), stored in Root objects, that specify the IntervalBox (a Cartesian product of intervals) and the “status” of a root. In this case, the package is able to guarantee, using an interval version of the Newton method, that there exists a root within each of those intervals, and that it is unique. Furthermore, any region of \mathbb{R}^2 that is not included in one of those intervals has been definitively excluded from the search.

Now let’s find the midpoints of the intervals and plot them:

midpoints = mid.([root.interval for root in rts])

xs = first.(midpoints)
ys = last.(midpoints)

using Plots; plotlyjs()

surface(-5:0.1:6, -6:0.1:6, (x,y)->f([x,y]))
scatter!(xs, ys, f.(midpoints))

The result is the following image:

35%20PM
35%20PM1114×994 188 KB

There are various usage examples in the README and in the examples.jl file. The documentation is currently out of date.

Please test the package with your difficult root-finding problems and open issues!

Currently the implementations are somewhat basic. There is lots of information in the literature about how to improve such root finders; this would make a great GSoC project. Please contact the first author if you are interested.

David P. Sanders
Luis Benet

41 Likes
2

Very nice! Can it also find roots that are connected, that is, roots that lie along a continuous curve in the xy-plane? I’m thinking specifically of finding the dispersion curves characterized by the Rayleigh-Lamb characteristic equations: https://en.wikipedia.org/wiki/Lamb_waves

2 Likes
3

Thanks!

Great question. No, for that you need IntervalConstraintProgramming.jl.

5 Likes
4

(And that package also needs work.)

5

Note that if you just want to plot implicit curves like that, without having any rigorous correctness guarantees, you can plot contours with the contour function in a Plots.jl.

1 Like
6

Do you mean that it needs work to solve those kinds of equations, or just that it needs to be updated/polished/refactored etc.?

7

How much work might be required to expand this to root finding on the complex plane?

1 Like
8

Great question – the answer is that it’s already implemented!

julia> using IntervalArithmetic, IntervalRootFinding

julia> x = -5..6
[-5, 6]

julia> Xc = Complex(x, x)
[-5, 6] + [-5, 6]*im

julia> f(z) = z^3 - 1
f (generic function with 1 method)

julia> rts = roots(f, Xc, Bisection)
7-element Array{IntervalRootFinding.Root{Complex{IntervalArithmetic.Interval{Float64}}},1}:
 Root([0.999511, 1.00019] + [-0.000183106, 0.000488282]*im, :unknown)
 Root([-0.500367, -0.499694] + [0.865905, 0.866578]*im, :unknown)
 Root([-0.500367, -0.499694] + [0.865234, 0.865906]*im, :unknown)
 Root([-0.501038, -0.500366] + [0.865905, 0.866578]*im, :unknown)
 Root([-0.499695, -0.499023] - [0.8656, 0.866272]*im, :unknown)
 Root([-0.500367, -0.499694] - [0.8656, 0.866272]*im, :unknown)
 Root([-0.500367, -0.499694] - [0.866271, 0.866944]*im, :unknown)

julia> rts = roots(f, rts, Newton)
3-element Array{IntervalRootFinding.Root{Complex{IntervalArithmetic.Interval{Float64}}},1}:
 Root([1, 1] + [-4.02199e-187, 5.17113e-187]*im, :unique)
 Root([-0.500001, -0.499999] + [0.866025, 0.866026]*im, :unique)
 Root([-0.500001, -0.499999] - [0.866025, 0.866026]*im, :unique)
4 Likes
9

I mean that it needs to be polished etc.
Feel free to try it out

The idea is to rewrite it using Cassette.jl, which could be another good GSoC project.

2 Likes
10

It is pretty naive – it unwraps the complex function as a function R^2->R^2 and uses the generic algorithm.

11

Functions of several complex variables are not yet implemented, but it should just be a question of extending the same technique.

12

Two questions:

(1) How do you guarantee that your sampling resolution is good enough?
(2) Is the package good enough to be considered a computer-assisted proof?

I am presuming no to both; because, naively you would need to take not just an evaluation oracle (including derivatives to arbitrary order via AD) but also non-local a-priori bounds. For example, if you told me that the gradient is Lipschitz with constant L over the domain, then I see how to use an approximate evaluation oracle + interval math + Newton to generate a computer-assisted proof. Without non-local a-priori bounds… well, for whatever evaluation points you have chosen, I can build another C^\infty function that matches at the evaluation points but has an extra zero somewhere else.

(You aren’t actually disassembling f, are you?)

Are there plans to improve the correctness up to computer-assisted-proof level? (presumably using Banach fixed-point instead of Brouwer, and just refusing to terminate in case of non-invertible jacobians at a zero)

That would be seriously cool for my non-computer-related research (but would need API changes, because the user must supply the a-priori bounds): There are many things that can be stated as root-finding where correctness/numerics is not an issue, but mathematical proof is (if I want to produce theorems, not pretty pictures).

13

Yes, the package should constitute a computer-assisted proof, since interval arithmetic, by definition, gives rigorous bounds for the image of a function f over a set X.

If 0 \notin f(X), where f(X) means “evaluate the function f using the rules of interval arithmetic, with input argument X” then the box X has been rigorously excluded from containing roots.
In general, interval arithmetic gives overestimates of the bounds.

e.g.

julia> using IntervalArithmetic

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

julia> X = -1..1
[-1, 1]

julia> f(X)
[-2, 3]

This is guaranteed to contain the true range, \{f(x): x \in X\}, even though in this case (and, in fact, in general, unless x only appears once in the expression of the function), the result is actually too large.

1 Like
14

If there is a non-invertible Jacobian, the package should probably return :unknown for the corresponding root, which mains that it is unable to guarantee the result.

15

Note that we do not use an evaluation oracle; these methods require that the function itself is evaluated completely using interval arithmetic. The code must be pure Julia and expressed as a combination of functions that are able to be evaluated using interval arithmetic.

e.g. for the moment you cannot use if...else; rather, you can reexpress the function using discontinuous functions such as abs, floor etc.

16

Ah, ok. Not that I would understand how interval math works internally; but small tests appear to confirm that the package is capable of finding very localized small, sharp spikes. Neat!

What about e.g.

f1(x)=2.5+sin(x) - 3.6 * exp(-200_000 * abs(x-1/pi))
roots(f1, 0..1, Bisection)
1-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([0.317382, 0.31836], :unknown)

This output is somwhat painful; the root-finder discovers that the box contains a zero and should be able to tell me that it contains even more zeros, or at least tell me the index of all the zeros in the box. Can I get this as output? If the solver manages to prove that there is a zero in some box then it presumably has computed the index, right?

17

There is an intro to interval arithmetic in the docs.

Note that in the nice, but difficult example you posted, the method has not confirmed that there is a zero; rather, it has proved that there is no zero outside of that interval.

The :unknown means that it can’t say for definite if there is one or there are more roots in that interval. The bisection method simply excludes boxes. If you want to prove that there is a root, you need to use (e.g.) the interval Newton method, which works for smooth functions (not sure what regularity is needed; probably at least C^1).

I tried replacing abs(x - 1/pi) by (x - 1/pi)^2 in your example, but even then Newton is unable to do anything, I guess because the spike is so narrow.

18

OK, I spoke too soon: you can just use BigFloat:

julia> f1(x) = 2.5 + sin(x) - 3.6 * exp(-200_000 * (x - 1/pi)^2);

julia> rts = roots(f1, big(-10..11), Newton, 1e-20)
2-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{BigFloat}},1}:
 Root([0.319419, 0.31942]₂₅₆, :unique)
 Root([0.317198, 0.317199]₂₅₆, :unique)

julia> rts
2-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{BigFloat}},1}:
 Root(Interval(3.194196639355450354064843277073749002133127335816266275112541721924471080255644e-01, 3.19419663935545073451199687688051252634626605118452346305176950426610552149445e-01), :unique)
 Root(Interval(3.171984202360467599768889278424171088018446102895213958259802499084915677675616e-01, 3.171984202360467854195005107499355294610187857277680252643125811521576337991869e-01), :unique)

julia> diam.([root.interval for root in rts])
2-element Array{BigFloat,1}:
 6.776263578034402712546580005437135696411132812500000000000000000000000000000000e-21
 6.776263578034402712546580005437135696411132812500000000000000000000000000000000e-21
19

(There seems to be some issue with the tolerance in Newton right now.)

20

By the way, I’m not sure what you mean by the “index”? (And no, it doesn’t calculate this, though if there’s some kind of formula for that it may be possible?)