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

**URL:** <https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515>\
**Category:** Community\
**Tags:** intervals, roots, scientific-computing\
**Created:** [March 5, 2018, 6:38pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515 "2018-03-05T18:38:42Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 5, 2018, 6:38pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/1 "2018-03-05T18:38:42Z")

</div>

We are excited to announce v0.2 of the [`IntervalRootFinding.jl`](https://github.com/JuliaIntervals/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](https://github.com/JuliaIntervals/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}:

```plaintext

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:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/1/2/129ed2d31dc4cdc51c43a45ffcf39f61da045708.jpeg)

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

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 5, 2018, 8:32pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/2 "2018-03-05T20:32:51Z")

</div>

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](https://en.wikipedia.org/wiki/Lamb_waves)

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 5, 2018, 8:38pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/3 "2018-03-05T20:38:24Z")

</div>

Thanks!

Great question. No, for that you need [`IntervalConstraintProgramming.jl`](https://github.com/dpsanders/IntervalConstraintProgramming.jl).

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 5, 2018, 8:47pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/4 "2018-03-05T20:47:57Z")

</div>

(And that package also needs work.)

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 5, 2018, 10:38pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/5 "2018-03-05T22:38:35Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 6, 2018, 11:31am UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/6 "2018-03-06T11:31:15Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![EP-Guy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ep-guy/32/3442_2.png) [@EP-Guy](https://discourse.julialang.org/u/EP-Guy)\
**Post date:** [March 7, 2018, 7:12pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/7 "2018-03-07T19:12:33Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 7, 2018, 7:47pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/8 "2018-03-07T19:47:26Z")

</div>

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

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 7, 2018, 7:50pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/9 "2018-03-07T19:50:42Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 7, 2018, 7:53pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/10 "2018-03-07T19:53:05Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 7, 2018, 7:57pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/11 "2018-03-07T19:57:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [March 7, 2018, 9:11pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/12 "2018-03-07T21:11:13Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 7, 2018, 9:21pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/13 "2018-03-07T21:21:10Z")

</div>

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
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.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 7, 2018, 9:22pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/14 "2018-03-07T21:22:00Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 7, 2018, 9:24pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/15 "2018-03-07T21:24:30Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [March 7, 2018, 10:23pm UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/16 "2018-03-07T22:23:14Z")

</div>

> [@dpsanders](#):
>
> Note that we do not use an evaluation oracle; these methods require that the function itself is evaluated completely using interval arithmetic.

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.

```julia
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?

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 8, 2018, 2:13am UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/17 "2018-03-08T02:13:40Z")

</div>

There is an intro to interval arithmetic in the [docs](https://juliaintervals.github.io/IntervalArithmetic.jl/stable/).

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.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 8, 2018, 2:16am UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/18 "2018-03-08T02:16:38Z")

</div>

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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 8, 2018, 2:17am UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/19 "2018-03-08T02:17:58Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 8, 2018, 5:18am UTC](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515/20 "2018-03-08T05:18:41Z")

</div>

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?)

[Next page](https://discourse.julialang.org/t/ann-intervalrootfinding-jl-for-finding-all-roots-of-a-multivariate-function/9515.md?page=2)