Determining the sign of a function over a given interval

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

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.

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

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

5 Likes

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

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

cc @mforets

What exactly do you need to do?

2 Likes

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!

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.

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.

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

2 Likes

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,

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:

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:

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

2 Likes

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

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.

2 Likes

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

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.

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

2 Likes

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

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]
2 Likes

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?

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

3 Likes

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?

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.

1 Like