# Determining the sign of a function over a given interval

**URL:** <https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231>\
**Category:** Optimization (Mathematical)\
**Created:** [January 11, 2020, 3:44am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231 "2020-01-11T03:44:38Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![ArjunNarayanan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arjunnarayanan/32/6475_2.png) [@ArjunNarayanan](https://discourse.julialang.org/u/ArjunNarayanan)\
**Post date:** [January 11, 2020, 3:44am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/1 "2020-01-11T03:44:38Z")

</div>

I am interested in determining the bounds of a function `f(x)` over a given interval `int`. Specifically, I want to determine if `f` is uniformly positive/negative over `int` or contains at least one root. `f` is continuous and finite so one of these three cases should occur.

Currently, I am using the excellent [TaylorModels.jl](https://github.com/JuliaIntervals/TaylorModels.jl) package by @dpsanders et al.

My approach is as follows:

- Compute a first order Taylor Model with bounded remainder of `f`.
- Evaluate the first order Taylor Model at the boundaries of `int` and add the remainder bound of the Taylor Model. Take the union of these two intervals.

Since the linear function is monotonic, the above interval bounds the range of `f` over `int`.

```julia
using TaylorModels

int = -1.0..1.0
x0 = mid(int)
f(x) = 3x^2 + 2x + 5
tm = TaylorModel1(1, x0, int)
ftm = f(tm)
left = ftm.pol(int.lo - x0) + ftm.rem
right = ftm.pol(int.hi - x0) + ftm.rem
bounds = union(left,right)
lower_bound = inf(bounds)
upper_bound = sup(bounds)

```

Now there are a few cases to consider:

- if `lower_bound > 0` then `f` is uniformly positive in `int`.
- if `upper_bound < 0` then `f` is uniformly negative in `int`.
- Otherwise, `f` possibly contains one (or more) roots in `int`.

To confirm that `f` has a zero crossing, I break the given interval into sub-intervals and repeat the above procedure. If I can determine at least one sub-interval where `f` is uniformly positive and at least one sub-interval where `f` is uniformly negative then there is guaranteed to be a zero crossing.

I was wondering if there is a more straightforward way of achieving this.

Specifically:

- Is there a smarter way of doing this for special functions? e.g. if `f(x)` is polynomial in `x`.
- Is there a way to use higher order Taylor Models? I am using first order because of monotonicity of the linear function.

I need the method to be general enough so as to work with multi-variables as well i.e. for functions `g(x,y)` and `h(x,y,z)` etc.

---

<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:** [January 11, 2020, 6:03am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/2 "2020-01-11T06:03:26Z")

</div>

> [@ArjunNarayanan](#):
>
> f(x) = 3x^2 + 2x + 5

The simplest thing to do is just use standard interval arithmetic, which happens to work in this case:

```julia
julia> using IntervalArithmetic

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

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

julia> f(X)
[3, 10]

```

To show that there is a zero crossing, if you know the function is continuous then it’s enough to find points `a` and `b` where `sign(f(a))` and `sign(f(b))` are opposite – points rather than intervals.

You can instead use the derivative (via automatic differentiation) and show that it is bounded away from 0:

```julia
julia> using ForwardDiff

julia> ForwardDiff.derivative(f, 2..3)
[14, 20]

```

This shows that the derivative is strictly positive over the interval, and hence the function is monotone there.

For more sophisticated root finding, check out `IntervalRootFinding.jl`:

```julia
julia> using IntervalRootFinding

julia> roots(f, X)
0-element Array{Root{Interval{Float64}},1}

```

This proves that there are no roots over the original interval.

cc @lbenet for Taylor models

---

<div class="post-metadata">

**Author:** ![ArjunNarayanan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arjunnarayanan/32/6475_2.png) [@ArjunNarayanan](https://discourse.julialang.org/u/ArjunNarayanan)\
**Post date:** [January 11, 2020, 6:42am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/3 "2020-01-11T06:42:48Z")

</div>

Thanks for your reply!

I initially started out with `IntervalArithmetic`, however the bounds with this approach were not tight. The input function `f(x)` can be a polynomial of much higher order than in my example – I think `IntervalArithmetic` would provide even looser bounds in this case (because of more arithmetic operations?)

My concern with the root finding approach is that it may not extend very well to a multivariate function `f(x,y)` etc.?

---

<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:** [January 11, 2020, 7:07am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/4 "2020-01-11T07:07:11Z")

</div>

Yes, `IntervalArithmetic` tends to produce rather loose bounds.  
IntervalRootFinding should extend to functions with a small number of variables.  
You can also try `IntervalOptimization`, which can find pretty good bounds on global maxima and minima.

There are various range bounding techniques in

[https://github.com/JuliaReach/RangeEnclosures.jl](https://github.com/JuliaReach/RangeEnclosures.jl)

cc @mforets

What exactly do you need to do?

---

<div class="post-metadata">

**Author:** ![ArjunNarayanan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arjunnarayanan/32/6475_2.png) [@ArjunNarayanan](https://discourse.julialang.org/u/ArjunNarayanan)\
**Post date:** [January 11, 2020, 8:25am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/5 "2020-01-11T08:25:40Z")

</div>

The exact problem statement is:

Given a multivariate polynomial p : \mathbb{R}^d \to \mathbb{R} and a hyper-rectangle R = (x\_1^L,x\_1^U) \times (x\_2^L, x\_2^U) \ldots (x\_d^L, x\_d^U), determine if p is uniformly positive, uniformly negative, or changes sign in R.

`RangeEnclosures.jl` looks very interesting, I will check it out, thanks!

---

<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:** [January 11, 2020, 9:19am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/6 "2020-01-11T09:19:20Z")

</div>

Right, that’s just the range bounding problem, or actually I guess the root-finding problem. It’s certainly a difficult problem in general, and gets of course exponentially harder as the dimension d increases.  
I think right now the most general solution is IntervalOptimisation.jl.

In the end you’re right that you need to split up the domain into pieces where you can prove what happens. That’s basically what all these packages do.

---

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [January 11, 2020, 7:22pm UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/7 "2020-01-11T19:22:15Z")

</div>

Hi @ArjunNarayanan, feel free to ask me if you need help with RangeEnclosures.jl. The pacakge lets you use different approaches available in the Julia ecosystem; on top of this, it will give an algorithm-agnostic branch-and-prune option through [GitHub - Kolaru/BranchAndPrune.jl: Branch and prune interface for Julia](https://github.com/Kolaru/BranchAndPrune.jl).

> Is there a smarter way of doing this for special functions? e.g. if `f(x)` is polynomial in `x` .

For a multivariate polynomial, you can try out polynomial optimization methods, available through the `:SumOfSquares` option (which runs [SumOfSquares.jl](https://github.com/JuliaOpt/SumOfSquares.jl)).

---

<div class="post-metadata">

**Author:** ![ArjunNarayanan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arjunnarayanan/32/6475_2.png) [@ArjunNarayanan](https://discourse.julialang.org/u/ArjunNarayanan)\
**Post date:** [January 12, 2020, 4:16am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/8 "2020-01-12T04:16:58Z")

</div>

Thanks for your reply @mforets.

I did take a look at `RangeEnclosures.jl`. I think combining `RangeEnclosures` with `BranchAndPrune` would solve my problem, but I’m not sure how to do it.

For example,

```julia
using RangeEnclosures
f(x) = (x - 0.2)*(x - 0.8)
int = 0..1
enclose(f,int) # [-0.640001, 0.160001]

```

In the above case, we cannot definitively say if `f` is uniformly positive, uniformly negative, or changes sign in `int`, and so we would have to:

> [@dpsanders](#):
>
> In the end you’re right that you need to split up the domain into pieces where you can prove what happens. That’s basically what all these packages do.

Is there a way to use `BranchAndPrune` in this process of subdivision?

---

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [January 12, 2020, 2:13pm UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/9 "2020-01-12T14:13:40Z")

</div>

> [@ArjunNarayanan](#):
>
> Is there a way to use `BranchAndPrune` in this process of subdivision?

For a start, here is a way to obtain a collection of intervals where `f` is monotonic using `BranchAndPrune`:

```julia
struct MonotoneSearch <: AbstractDepthFirstSearch{Interval{Float64}}
    f::Function
    initial::Interval
    tol::Float64
    algorithm::Symbol
end

is_monotonic(x::Interval) = inf(x)*sup(x) > 0

function BranchAndPrune.process(search::MonotoneSearch, interval)
    y = enclose(search.f, interval, search.algorithm)

    if is_monotonic(y) # sign doesn't change
        return :store, interval
    elseif diam(y) < search.tol
        return :discard, interval
    else
        return :bisect, interval
    end
end

BranchAndPrune.bisect(::MonotoneSearch, interval) = bisect(interval)

function run_search(f, interval; tol=1/256, alg=:IntervalArithmetic)
    search = MonotoneSearch(f, interval, tol, alg)

    local endtree = nothing

    for working_tree in search
        endtree = working_tree
    end

    return endtree
end

```

```julia
julia> tree = run_search(f, int)
Working tree with 29 elements of type Interval{Float64}
Indices: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 20, 21, 22, 23, 24, 25, 26, 27, 28, 30]
Structure:
  [1] Node with children [2, 3]
    [2] Node with children [17, 18]
      [17] Node with children [19, 20]
        [19] Leaf (:final) with data [0, 0.122094]
        [20] Node with children [21, 22]
          [21] Leaf (:final) with data [0.122093, 0.183617]
          [22] Node with children [23, 24]
            [23] Node with children [25, 26]
              [25] Leaf (:final) with data [0.183616, 0.198997]
              [26] Node with children [27, 28]
                [27] Node with children [30]
                  [30] Leaf (:final) with data [0.202841, 0.206747]
                [28] Leaf (:final) with data [0.206746, 0.214619]
            [24] Leaf (:final) with data [0.214618, 0.24611]
      [18] Leaf (:final) with data [0.246109, 0.496094]
    [3] Node with children [4, 5]
      [4] Leaf (:final) with data [0.496093, 0.746079]
      [5] Node with children [6, 7]
        [6] Node with children [8, 9]
          [8] Node with children [10, 11]
            [10] Leaf (:final) with data [0.746078, 0.777081]
            [11] Node with children [12, 13]
              [12] Leaf (:final) with data [0.77708, 0.792703]
              [13] Node with children [14, 15]
                [14] Node with children [16]
                  [16] Leaf (:final) with data [0.792702, 0.796609]
                [15] Leaf (:final) with data [0.800574, 0.808571]
          [9] Leaf (:final) with data [0.80857, 0.872048]
        [7] Leaf (:final) with data [0.872047, 1]

# sanity check
julia> all(is_monotonic(enclose(f, x)) for x in data(tree))
true

```

The multivariate case can be treated similarly, see this issue: [Add branch and prune method · Issue #6 · JuliaReach/RangeEnclosures.jl · GitHub](https://github.com/JuliaReach/RangeEnclosures.jl/issues/6)

---

<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:** [January 13, 2020, 7:09am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/10 "2020-01-13T07:09:32Z")

</div>

Here’s an alternative solution. We can change the question to “find the region where the function is positive”:

```julia
julia> using IntervalArithmetic, IntervalConstraintProgramming, ModelingToolkit

julia> f(x) = (x - 0.2)*(x - 0.8)
f (generic function with 1 method)

julia> vars = @variables x
(x,)

julia> S = Separator(x, f(x) > 0);

julia> X = IntervalBox(-5..5)
[-5, 5]

julia> p = pave(S, X, 0.1)
Paving:
- tolerance ϵ = 0.1
- inner approx. of length 5
- boundary approx. of length 2

julia> p.inner
5-element Array{IntervalBox{1,Float64},1}:
 [2.46078, 5]
 [1.20109, 2.46079]
 [0.8, 1.2011]
 [-0.0390625, 0.200001]
 [-5, -0.0390625]

julia> p.boundary
2-element Array{IntervalBox{1,Float64},1}:
 [0.8, 0.800001]
 [0.2, 0.200001]

```

This uses a branch and bound technique together with interval constraint propagation.

`p.inner` are intervals that are proved to satisfy the inequality. `p.boundary` are intervals where it can’t prove if they satisfy the inequality or not. Intervals that are proved _not_ to satisfy the inequality have been removed.

This is a very powerful technique, but starts to break down when the dimension is too high if the sets have complicated shapes, since it represents complicated shapes as unions of boxes.

---

<div class="post-metadata">

**Author:** ![ArjunNarayanan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arjunnarayanan/32/6475_2.png) [@ArjunNarayanan](https://discourse.julialang.org/u/ArjunNarayanan)\
**Post date:** [January 16, 2020, 6:06am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/11 "2020-01-16T06:06:03Z")

</div>

Just wanted to build on @mforets answer specifically for the original question in case someone is interested in a similar problem in the future,

```julia
using BranchAndPrune, RangeEnclosures, IntervalArithmetic
using BenchmarkTools, Plots

# Define the search type as mutable so we can change the value of
# flags in order to terminate the search.
# Breadth first search was more applicable in this search problem.
mutable struct BoundSearch{T} <: AbstractBreadthFirstSearch{Interval{T}}
    f::Function
    initial::Interval
    algorithm::Symbol
    tol::Real
    found_positive::Bool
    found_negative::Bool
    breached_tolerance::Bool
    function BoundSearch{T}(f::Function, initial::Interval, algorithm::Symbol, tol::Real) where {T}

        if tol < 0
            throw(ArgumentError("Tolerance must be a positive number"))
        end
        new{T}(f,initial,algorithm,tol,false,false,false)
    end
end

function BranchAndPrune.process(search::BoundSearch, interval::Interval)

    f_range = enclose(search.f, interval, search.algorithm)

    if (search.found_positive && search.found_negative) || (search.breached_tolerance)
        # No further processing required in this case
        return :discard, interval
    elseif inf(f_range) > 0.0
        # We have found a sub-interval where the function is
        # uniformly positive
        search.found_positive = true
        return :store, interval
    elseif sup(f_range) < 0.0
        # We have found a sub-interval where the function is
        # uniformly negative
        search.found_negative = true
        return :store, interval
    elseif diam(interval) < search.tol
        search.breached_tolerance = true
        return :discard, interval
    else
        return :bisect, interval
    end
end

BranchAndPrune.bisect(::BoundSearch, interval) = bisect(interval,0.5)

function run_search(f, interval, algorithm, tol)

    search = BoundSearch{Float64}(f, interval, algorithm, tol)
    local endtree = nothing
    for working_tree in search
        endtree = working_tree
    end
    return endtree, search
end

"""
    sign(f::Function, int::Interval)
return
- `+1` if `f` is uniformly positive on `int`
- `-1` if `f` is uniformly negative on `int`
- `0` if `f` has at least one zero crossing in `int` (f assumed continuous)
"""
function Base.sign(f::Function, int::Interval; algorithm = :IntervalArithmetic, tol = 1e-2)
    tree, search = run_search(f,int,algorithm,tol)
    if search.found_positive && search.found_negative
        return 0
    elseif search.breached_tolerance
        error("Search tolerance not tight enough")
    elseif search.found_positive
        return 1
    else search.found_negative
        return -1
    end
end

```

The major change is to declare our search type as mutable so we can set the value of certain flags based on the search results – this allows us to terminate the search once certain criteria are fulfilled.  
Further, we subtype `AbstractBreadthFirstSearch`: the status of an interval around a root of `f` is undetermined because the bounds will always be `[+ve number, -ve number]`. Thus, if there is a root, depth first search will always exceed the given tolerance which is undesirable.

```julia
f(x) = (x - 0.0)*(x - 0.1)*(x - 1.0)
int = 0 .. 1
s = @btime sign(f,int) # 32.878 μs (415 allocations: 27.25 KiB)

```

![example_poly](https://global.discourse-cdn.com/julialang/original/3X/b/a/ba532b1220ef33e92f379e33c36e0d097f730d8e.png)

---

<div class="post-metadata">

**Author:** ![lbenet](https://avatars.discourse-cdn.com/v4/letter/l/35a633/32.png) [@lbenet](https://discourse.julialang.org/u/lbenet)\
**Post date:** [February 13, 2020, 10:20pm UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/12 "2020-02-13T22:20:09Z")

</div>

My two cents contribution on this, using TaylorModels.jl:

```julia
julia> using TaylorModels

julia> f(x) = (x - 0.0)*(x - 0.1)*(x - 1.0)
f (generic function with 1 method)

julia> Dom = 0..1 # Domain of interest
[0, 1]

julia> xTM = TaylorModel1( 3, 0..0, Dom) # independent variable or order/degree 3
 [1, 1] t + [0, 0]

julia> y = f(xTM) # Taylor model of `f` on the domain `Dom`
 [0.1, 0.100001] t + [-1.10001, -1.09999] t² + [1, 1] t³ + [0, 0]

julia> evaluate(y, Dom) # range of `f` over `Dom`; same result as a naive Interval evaluation
[-1.00001, 0.100001]

julia> ivs = mince(Dom, 32); # divide the domain in 32 intervals

julia> hull(evaluate.(y, ivs)) # hull of the ranges of each small interval
[-0.141224, 0.00422364]

```

---

<div class="post-metadata">

**Author:** ![ArjunNarayanan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arjunnarayanan/32/6475_2.png) [@ArjunNarayanan](https://discourse.julialang.org/u/ArjunNarayanan)\
**Post date:** [February 14, 2020, 2:02am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/13 "2020-02-14T02:02:25Z")

</div>

This is very useful. While I think `RangeEnclosures.jl` is a great package, I think it is a little more than what I need for my purposes. There are also a great many dependencies that I would rather not want to introduce into my own package. A pure `TaylorModels` solution is very useful.

A follow up question: In the case that `f(x)` is a polynomial (as it is in my/your example), is there a faster method within the `TaylorModels` package to bound the range? I have briefly gone through the source code and come across a few remainder bounding methods and was wondering if I could directly apply any one of those to my function?

---

<div class="post-metadata">

**Author:** ![lbenet](https://avatars.discourse-cdn.com/v4/letter/l/35a633/32.png) [@lbenet](https://discourse.julialang.org/u/lbenet)\
**Post date:** [February 14, 2020, 3:33am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/14 "2020-02-14T03:33:44Z")

</div>

If `f(x)` is a polynomial, mincing is a safe way to proceed. As I showed above, using a `TaylorModel1` of order 3 gives the exact representation (zero remainder), but naive evaluation overestimates the range, essentially in the same way pure interval arithmetic does. Yet, direct mincing converges slowly.

For the example above (a one dimensional polynomial), you can also do the following:

```julia
julia> TaylorModels.bound_taylor1(y, Dom)
[-0.126229, 0.00237665]

```

In this case, the range is much tighter. The method uses `IntervalRootFinding.jl` under the hood to approximate the roots of the derivative (maxima and minima) and then get tighter bounds. For polynomials _in one variable_ I think the method is stable; we haven’t extended it for several variables.

---

<div class="post-metadata">

**Author:** ![ArjunNarayanan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arjunnarayanan/32/6475_2.png) [@ArjunNarayanan](https://discourse.julialang.org/u/ArjunNarayanan)\
**Post date:** [February 15, 2020, 12:08am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/15 "2020-02-15T00:08:06Z")

</div>

Thanks! I am definitely also interested in the multivariate case.  
I think getting the (dramatically) tighter bound that you have obtained would be quite useful for the application I am working on; mainly from an efficiency perspective. As you point out, a mincing strategy does take some time to converge.  
Would the multivariate case be worth a PR from my side? Do you have a strategy in mind that I could potentially look into?

---

<div class="post-metadata">

**Author:** ![lbenet](https://avatars.discourse-cdn.com/v4/letter/l/35a633/32.png) [@lbenet](https://discourse.julialang.org/u/lbenet)\
**Post date:** [February 15, 2020, 1:01am UTC](https://discourse.julialang.org/t/determining-the-sign-of-a-function-over-a-given-interval/33231/16 "2020-02-15T01:01:25Z")

</div>

> [@ArjunNarayanan](#):
>
> Would the multivariate case be worth a PR from my side? Do you have a strategy in mind that I could potentially look into?

Yes, a PR is certainly welcome! Regarding a strategy, I would first try finding the maxima and minima using IntervalRootFinding.jl, or perhaps a combination with BranchAndPrune.jl.