# How can I check if a function is always positive?

**URL:** <https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405>\
**Category:** General Usage\
**Tags:** question\
**Created:** [October 15, 2020, 2:26am UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405 "2020-10-15T02:26:50Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)\
**Post date:** [October 15, 2020, 2:26am UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/1 "2020-10-15T02:26:50Z")

</div>

Suppose I have a user-defined mathematical function, say

```julia
f(x::Real) = x*exp(x)

```

I would like to check if f(x) \> 0 for all x \in [x\_l, x\_u] \subset \mathbb{R}.

The best ideas I have are either

```julia
x = xₗ:0.01:xᵤ
if any(f.(x) .≤ 0)
    ErrorException("Input function must be positive for x ∈ [$xₗ, $xᵤ].")
end 

```

or this option (which only checks dynamically, when I want it checked statically (if possible))

```julia
function f_internal(x)
    val = f(x)
    if val ≤ 0
        ErrorException("Input function must be positive for x ∈ [$xₗ, $xᵤ].") |> throw
    else
        return val
    end
end

```

but is there a better/Julian way to do this?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 15, 2020, 2:34am UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/2 "2020-10-15T02:34:09Z")

</div>

You should check out IntervalArithmatic

---

<div class="post-metadata">

**Author:** ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)\
**Post date:** [October 15, 2020, 2:42am UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/3 "2020-10-15T02:42:03Z")

</div>

Oh that’s perfect. Thanks!

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 15, 2020, 2:51am UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/4 "2020-10-15T02:51:44Z")

</div>

The one warning is that the basic Interval methods can be slow if your function is complicated/ is badly behaved in some technical ways. if IntervalArithmatic doesn’t work well, you might want to try TaylorModels which is a more sophisticated (but higher overhead) technique.

---

<div class="post-metadata">

**Author:** ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)\
**Post date:** [October 15, 2020, 2:52am UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/5 "2020-10-15T02:52:40Z")

</div>

Oh thanks! If I have issues with IntervalArithmetic, then yes I’ll check out TaylorModels.

---

<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:** [October 15, 2020, 4:03am UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/6 "2020-10-15T04:03:58Z")

</div>

Here’s an example using [GitHub - JuliaIntervals/IntervalArithmetic.jl: Rigorous floating-point calculations using interval arithmetic in Julia](https://github.com/JuliaIntervals/IntervalArithmetic.jl).

The main method / goal of interval arithmetic is to calculate the **range** of a function f over an input set X, i.e. the set \mathcal{R}(f; X) := \{f(x): x \in X \}. Calculating this exactly is hard (impossible in general, equivalent to global optimisation of the function), but interval arithmetic gives a cheap way of calculating an **enclosure** Y of this range, i.e. an interval that is guaranteed to contain it: \mathcal{R}(f; X) \subseteq Y.

The idea is to split up a complicated function f into simple pieces and make sure you bound the range for each piece. Putting the pieces together gives you an enclosure for the complicated function (the “fundamental theorem of interval arithmetic”), but which is in general an _over_-estimate.

```julia
julia> using IntervalArithmetic

julia> f(x) = x * exp(x)
f (generic function with 1 method)

julia> X = interval(1, 2)
[1, 2]

julia> f(X)
[2.71828, 14.7782]

julia> f(X).lo > 0
true

```

This constitutes a _rigorous proof_ that f(x) \> 0 for all x \in X, where X = [1, 2] is the interval of all real numbers between 1 and 2.

In fact you can also do things like

```julia
julia> f(0..∞)
[0, ∞]

```

to prove this for _any_ x \ge 0!

But this technique has “some” (many) limitations, where it can’t prove things that you know to be true, for example

```julia
julia> g(x) = x^2 - 2x

julia> g(100..∞)
[-∞, ∞]

```

This just says that the true range of the function over that input set is _contained in_ (i.e. is a subset of) the given result.

In this case we (humans) actually know that x^2 is much bigger than 2x in that interval, and hence the function is positive there. In this particular case you could deduce that by calculating the derivative, e.g. using `ForwardDiff.jl` (which works out of the box with intervals!). But in general I believe it is a hard problem. (But I am _very willing_ to be proved wrong here!)

---

<div class="post-metadata">

**Author:** ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)\
**Post date:** [October 15, 2020, 4:13am UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/7 "2020-10-15T04:13:45Z")

</div>

The concept is awesome, and I’m loving the implementation!

The example with g(x) = x^2 - 2x scares me though, because g is a composition of elementary functions and `IntervalArithmetic` doesn’t process that it’s positive for all x \in [100, \infty).

So, I’ll try and be careful with the implementation I have in my code. Luckily, I just need this functionality for checking output range rather than complicated and involved transformations on the input range.

Thanks @dpsanders!

---

<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:** [October 15, 2020, 4:20am UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/8 "2020-10-15T04:20:40Z")

</div>

In general the over-estimation of the range comes when a variable is repeated more than once in an expression, as is the case for x^2 - 2x. This is due to the definition of subtraction of two intervals as

X - Y := \{x - y: x \in X, y \in Y \}

So that if e.g. X = [0, 1] then X - X = [-1, 1] instead of 0. This is known as the **dependency effect**. It is difficult to get rid of, unfortunately.

In the case of g(x) = x^2 - 2x you can do the following:

```julia
julia> g(X.lo)
9800.0

julia> ForwardDiff.derivative(g, X)
[198, ∞]

```

Since the derivative is positive everywhere in X, the function is increasing. Since it starts out positive, it remains positive.

But then

```julia
julia> h(x) = x^3 - 3x^2
h (generic function with 1 method)

julia> ForwardDiff.derivative(h, X)
[-∞, ∞]

```

so in that case you would need to go to the _second_ derivative.

For a complicated transcendental function I’m not sure what you would do.

---

<div class="post-metadata">

**Author:** ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)\
**Post date:** [October 15, 2020, 4:23am UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/9 "2020-10-15T04:23:28Z")

</div>

Ah, that makes sense. Utilize the derivatives to enable checking of the range, in a sense, semi-manually constructing the range of the function in question.

@dpsanders or anyone, how well does `IntervalArithmetic` play with interpolated functions such as solutions to `DifferentialEquations`?

---

<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:** [October 15, 2020, 4:29am UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/10 "2020-10-15T04:29:54Z")

</div>

Good question. I’m guessing they’re just polynomials and you are probably OK?  
But the output of a solver from `DifferentialEquations` is never rigorous.

If you want guaranteed correct (rigorous) bounds for solutions of ODEs you need `TaylorModels.jl`.

---

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [October 15, 2020, 10:03am UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/11 "2020-10-15T10:03:40Z")

</div>

Btw. @dpsanders made a really nice hands-on tutorial on this topic: [JuliaCon 2020 | Calculating with Sets: Interval Methods in Julia - YouTube](https://www.youtube.com/watch?v=LAuRCy9jUU8)

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [October 16, 2020, 11:22am UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/12 "2020-10-16T11:22:56Z")

</div>

It also matters how you write it:

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

julia> f1(100..Inf)
[-∞, ∞]

julia> f2(x) = x*(x - 2) # mathematically identical to f1
f2 (generic function with 1 method)

julia> f2(100..Inf)
[9800, ∞]

```

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [October 16, 2020, 12:59pm UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/13 "2020-10-16T12:59:01Z")

</div>

Why not use BlackBoxOptim.jl to o find the global minimum, and compare with `0`?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 16, 2020, 1:11pm UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/14 "2020-10-16T13:11:41Z")

</div>

> [@kapple](#):
>
> Suppose I have a user-defined mathematical function […] I would like to check if f(x) \> 0 for all x \in [x\_l, x\_u] \subset \mathbb{R}.

One alternative here is to use [ApproxFun.jl](https://github.com/JuliaApproximation/ApproxFun.jl), e.g.

```julia-auto
julia> using ApproxFun

julia> checkpos(f, xₗ, xᵤ) = minimum(Fun(f, xₗ..xᵤ)) > 0

julia> checkpos(x -> x*exp(x), 3, 4)
true

julia> checkpos(x -> x*exp(x), -2, 2)
false

julia> checkpos(x -> x^2 - 2x, 100, 200)
true

```

This works by constructing a polynomial approximant to the function, to nearly machine precision, and then finding the minimum of that polynomial.

Unlike interval arithmetic, it doesn’t have a dependency problem where it drastically overestimates the range. On the other hand, it might give a false positive/negative if roundoff errors cause the function approximation to slightly cross the axis, and it only allows you to check finite (and reasonably small) intervals.

---

<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:** [October 16, 2020, 3:15pm UTC](https://discourse.julialang.org/t/how-can-i-check-if-a-function-is-always-positive/48405/15 "2020-10-16T15:15:35Z")

</div>

That’s a good point, you can also use IntervalOptimisation.jl and IntervalConstraintProgramming.jl to do this, at least over finite intervals.