# Interval arithmetic - computation time

**URL:** <https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633>\
**Category:** Numerics\
**Tags:** intervals\
**Created:** [September 6, 2018, 5:58pm UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633 "2018-09-06T17:58:07Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![OlivierHnt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olivierhnt/32/6227_2.png) [@OlivierHnt](https://discourse.julialang.org/u/OlivierHnt)\
**Post date:** [September 6, 2018, 5:58pm UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/1 "2018-09-06T17:58:08Z")

</div>

Hi, I just started using the IntervalArithmetic package for Julia. For the same computations, Matlab (using IntLab) performs extremely well in comparison with Julia.

The implementation has ben adapted, that is the computations are vectorized in Matlab but not in Julia. Thus, as expected if I compare the code with no interval arithmetic, Julia is faster than Matlab.

Is there a specific and more efficient way of computing when using IntervalArithmetic ? Currently, my code without interval arithmetic and with interval arithmetic are the same except that I ensure that all the variables are Interval{Float64}.

Thank you for your help!

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [September 6, 2018, 6:06pm UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/2 "2018-09-06T18:06:30Z")

</div>

It’s really hard to say anything at all without some idea of what your code is.

Do you have an example of a (relatively short) function which exhibits slow performance with intervals in Julia?

---

<div class="post-metadata">

**Author:** ![OlivierHnt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olivierhnt/32/6227_2.png) [@OlivierHnt](https://discourse.julialang.org/u/OlivierHnt)\
**Post date:** [September 6, 2018, 6:28pm UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/3 "2018-09-06T18:28:09Z")

</div>

Your reply then suggests that there is no function a priori costly that should be avoided. Good to know!

About the code, here is one:

```julia
using LinearAlgebra
using IntervalArithmetic
using BenchmarkTools
using Base.Threads

function fun1(x,n,l)
    sqrt2 = sqrt(2)
    F = Matrix{Float64}(I, n, n)
    @inbounds for j in 1:n
        @inbounds @threads for i in 1:n
            @inbounds for k in 0:l-1
                theta = k*2*pi/l
                costheta = cos(theta)
                cos2theta = cos(2*theta)
                xisquare, xjsquare, xixj = x[i]^2, x[j]^2, x[i]*x[j]
                if i!=j
                    d = (xisquare + xjsquare - 2*xixj*costheta)^(1/2)
                    F[i,i] += -(4*xisquare + xjsquare + xjsquare*cos2theta)/d^5
                    F[i,j] += 4*(xisquare + xjsquare)*costheta/d^5
                elseif i==j && k>0
                    F[i,j] += -1/x[i]^3/(1-costheta)^(1/2)
                end
            end
        end
    end
    return F
end

function fun2(ix,n,l)
    sqrt2 = @interval(sqrt(2)) # this should be safe
    iF = Matrix{Interval{Float64}}(I, n, n)
    @inbounds for j in 1:n
        @inbounds @threads for i in 1:n
            @inbounds for k in 0:l-1
                itheta = k*2*@interval(pi)/l
                icostheta = cos(itheta)
                icos2theta = cos(2*itheta)
                ixisquare, ixjsquare, ixixj = ix[i]^2, ix[j]^2, ix[i]*ix[j]
                if i!=j
                    id = (ixisquare + ixjsquare - 2*ixixj*icostheta)^(1/2)
                    iF[i,i] += -(4*ixisquare + ixjsquare + ixjsquare*icos2theta)/id^5
                    iF[i,j] += (ixisquare + ixjsquare)*icostheta/id^5
                elseif i==j && k>0
                    iF[i,j] += -1/ix[i]^3/(1-icostheta)^(1/2)
                end
            end
        end
    end
    return iF
end

n=20
l=20

x = float(reshape(collect(1:n),n,1))
ix = convert(Array{Interval{Float64},2}, x)

@btime fun1(x,n,l)

@btime fun2(ix,n,l)

```

The code on which I noticed the time difference is a more elaborated version of this one but the loops are the same (I just remove a lot of terms to make it easier to read).  
Using IntervalArithmetic, the runtime is almost multiplied by a factor 200.

---

<div class="post-metadata">

**Author:** ![ExpandingMan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/expandingman/32/866_2.png) [@ExpandingMan](https://discourse.julialang.org/u/ExpandingMan)\
**Post date:** [September 6, 2018, 6:30pm UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/4 "2018-09-06T18:30:58Z")

</div>

> [@OlivierHnt](#):
>
> Your reply then suggests that there is no function a priori costly that should be avoided.

At the risk of putting words into the core devs mouths, I believe that is one of the core tenants of the design philosophy of Julia (this may or may not apply to IntervalArithmetic).

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [September 6, 2018, 6:55pm UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/5 "2018-09-06T18:55:28Z")

</div>

Hm, it looks like there is a performance issue with literal powers of intervals. For example:

```julia
julia> x = Interval(1.0, 1.0)
[1, 1]

julia> @btime $x * $x
  8.741 ns (0 allocations: 0 bytes)
[1, 1]

julia> @btime $x^2
  1.393 μs (51 allocations: 2.03 KiB)
[1, 1]

julia> @btime sqrt($x)
  56.531 ns (0 allocations: 0 bytes)
[1, 1]

julia> @btime $x^(1/2)
  1.940 μs (67 allocations: 2.78 KiB)
[1, 1]

```

It looks like the relevant issue, with some workarounds, is here: [Power ^ allocates · Issue #103 · JuliaIntervals/IntervalArithmetic.jl · GitHub](https://github.com/JuliaIntervals/IntervalArithmetic.jl/issues/103)

Just by removing all the literal powers from your sample code, I was able to get a speedup of about 10x.

---

<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:** [September 6, 2018, 7:07pm UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/6 "2018-09-06T19:07:52Z")

</div>

> [@rdeits](#):
>
> Just by removing all the literal powers from your sample code, I was able to get a speedup of about 10x.

`sqrt(x)` is almost always going to be faster than `x^0.5` anyway, because the latter dispatches to generic code whereas the former is optimized for square roots.

There are probably a few other small optimizations, e.g. `-1/ix[i]^3/(1-icostheta)^(1/2)` could be replaced by `inv(ix[i]*ix[i]*ix[i] * sqrt((1-icostheta))`, to only do the division once (and with specialized `inv` code rather than a generic division of `1`).

---

<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:** [September 6, 2018, 8:04pm UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/7 "2018-09-06T20:04:50Z")

</div>

Hi,

As it has been pointed out, there are some optimizations related to `sqrt`, avoid some repeated calculations some constants that you can move out from the loops and divisions, and powers are an issue.

```julia
function fun3(ix,n,l)
    sqrt2 = @interval(sqrt(2)) # this should be safe
    piI = @interval(pi)
    iF = Matrix{Interval{Float64}}(I, n, n)
    @inbounds for j in 1:n
       @inbounds @threads for i in 1:n
           @inbounds for k in 0:l-1
               itheta = k*2*piI/l
               icostheta = cos(itheta)
               icos2theta = cos(2*itheta)
               ixisquare, ixjsquare, ixixj = ix[i]^2, ix[j]^2, ix[i]*ix[j]
               if i!=j
                   id = sqrt(ixisquare + ixjsquare - 2*ixixj*icostheta)
                   id5 = id^5
                   iF[i,i] += -(4*ixisquare + ixjsquare + ixjsquare*icos2theta)/id5
                   iF[i,j] += (ixisquare + ixjsquare)*icostheta/id5
               elseif i==j && k>0
                   iF[i,j] += -1/(ix[i]^3 * sqrt(1-icostheta))
               end
           end
       end
    end
    return iF
end

```

In my machine with Julia 0.7, I get the current benchmarks:

```julia
@btime fun1(x,n,l)
  5.363 ms (238781 allocations: 3.65 MiB)

@btime fun2(ix,n,l)
  179.564 ms (2656832 allocations: 96.55 MiB)

@btime fun3(ix,n,l)
  104.914 ms (1632388 allocations: 57.50 MiB)

```

---

<div class="post-metadata">

**Author:** ![tshort](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tshort/32/43_2.png) [@tshort](https://discourse.julialang.org/u/tshort)\
**Post date:** [September 6, 2018, 8:09pm UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/8 "2018-09-06T20:09:13Z")

</div>

`fun3` didn’t fix all the powers (the biggest issue). I think you can replace `ix[i]^2` by `pow(ix[i], 2)` or write out the multiplication.

---

<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:** [September 6, 2018, 8:09pm UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/9 "2018-09-06T20:09:23Z")

</div>

And using `pow(x, n)` instead of `x^n`, you get better timings:

```julia
julia> function fun4(ix,n,l)
    sqrt2 = @interval(sqrt(2)) # this should be safe
    piI = @interval(pi)
    iF = Matrix{Interval{Float64}}(I, n, n)
    @inbounds for j in 1:n
       @inbounds @threads for i in 1:n
           @inbounds for k in 0:l-1
               itheta = k*2*piI/l
               icostheta = cos(itheta)
               icos2theta = cos(2*itheta)
               ixisquare, ixjsquare, ixixj = pow(ix[i],2), pow(ix[j],2), ix[i]*ix[j]
               if i!=j
                   id = sqrt(ixisquare + ixjsquare - 2*ixixj*icostheta)
                   id5 = pow(id,5)
                   iF[i,i] += -(4*ixisquare + ixjsquare + ixjsquare*icos2theta)/id5
                   iF[i,j] += (ixisquare + ixjsquare)*icostheta/id5
               elseif i==j && k>0
                   iF[i,j] += -1/(pow(ix[i],3) * sqrt(1-icostheta))
               end
           end
       end
    end
    return iF
end

julia> @btime fun4(ix,n,l)
  17.375 ms (231206 allocations: 7.06 MiB)

```

---

<div class="post-metadata">

**Author:** ![OlivierHnt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olivierhnt/32/6227_2.png) [@OlivierHnt](https://discourse.julialang.org/u/OlivierHnt)\
**Post date:** [September 6, 2018, 8:30pm UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/10 "2018-09-06T20:30:23Z")

</div>

Thank you all for your suggestions!  
I have been able to correct my code with great results!

---

<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:** [September 6, 2018, 9:09pm UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/11 "2018-09-06T21:09:14Z")

</div>

Can you post your benchmarks to compare with IntLab?

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [September 6, 2018, 10:08pm UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/12 "2018-09-06T22:08:58Z")

</div>

Also, if you see performance issues with `@threads` (for example, if removing it actually makes your code faster), then you’re probably running into [performance of captured variables in closures · Issue #15276 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/15276) . There’s some discussion of the issue here: [Parallelizing for loop in the computation of a gradient - #7 by tkoolen](https://discourse.julialang.org/t/parallelizing-for-loop-in-the-computation-of-a-gradient/9154/7) as well.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [September 7, 2018, 3:13am UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/13 "2018-09-07T03:13:23Z")

</div>

FWIW only the outermost `@inbounds` is needed, it applies to the whole block, including the inner loops.

---

<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:** [September 7, 2018, 8:31am UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/14 "2018-09-07T08:31:10Z")

</div>

I’m a bit confused about what your code is doing here. What is the purpose of

> [@OlivierHnt](#):
>
> ```julia
> sqrt2 = @interval(sqrt(2))
> 
> ```

?  
It’s not used. And why is this here

> [@OlivierHnt](#):
>
> ```julia
> itheta = k*2*@interval(pi)/l
> 
> ```

?

Why do you write two different functions for ordinary floats and intervals instead of a single one that handles both?

I would have expected that

```julia
all(fun1(x, n, l) .∈ fun2(ix, n, l))

```

to be true, but it’s not.

If I write a single function, like this, I get something more understandable:

```julia
mypow(x, n::Integer) = x^n
mypow(x::Interval, n::Integer) = pow(x, n)

function fun0((x::AbstractVector, l::Integer)
    F = diagm(0 => x)
    @inbounds for j in eachindex(x) # dangerous to use 1:n where n is an input, since you have @inbounds
        xj2 = x[j] * x[j]
        Threads.@threads for i in eachindex(x)
            xi2 = x[i] * x[i]
            xixj = x[i] * x[j]
            for k in 0:l-1
                θ = 2π * k / l # I did not make this an interval. Why should it be an interval?
                cosθ = cos(θ)
                cos2θ = cos(2θ)
                if i != j
                    d = 1 / mypow(sqrt(xi2 + xj2 - 2xixj * cosθ), 5)
                    F[i,i] += -(4xi2 + xj2 + xj2 * cos2θ) * d
                    F[i,j] += 4 * (xi2 + xj2) * cosθ * d
                elseif k > 0
                    F[i, j] -= 1 / (sqrt(1 - cosθ) * mypow(x[i], 3))
                end
            end
        end
    end
    return F
end

```

Then call it with:

```julia
julia> x = float.(1:n); # make vector not matrix!

julia> xi = interval.(x); # much simpler notation (also, use vector)

julia> all(fun0(x, 20) .∈ fun0(xi, 20)) # this is what I would have expected
true

julia> @btime fun0(x, 20);
  4.941 ms (68758 allocations: 1.03 MiB)

julia> @btime fun0(xi, 20); # pretty fast
  6.383 ms (80461 allocations: 2.33 MiB)

```

**Edit:** Forgot to interpolate `x` and `xi` into `@btime`, but it makes no difference.  
**Edit2:** For comparison with the old implementation for non-intervals:

```julia
julia> @btime fun1($x, 20, 20);
  6.335 ms (80782 allocations: 1.22 MiB)

```

---

<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:** [September 7, 2018, 8:35am UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/15 "2018-09-07T08:35:43Z")

</div>

BTW, since the matrices are symmetrical, you can make it faster by exploiting that.

---

<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:** [September 7, 2018, 9:14am UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/16 "2018-09-07T09:14:31Z")

</div>

Exploiting symmetry:

```julia
function funS(x::AbstractVector, l::Integer)
    F = diagm(0 => x)
    @inbounds for j in eachindex(x)
        xj2 = x[j] * x[j]
        Threads.@threads for i in 1:j
            xi2 = x[i] * x[i]
            xixj = x[i] * x[j]
            for k in 0:l-1
                θ = 2π * k / l
                cosθ = cos(θ)
                cos2θ = cos(2θ)
                if i != j
                    d = 1 / mypow(sqrt(xi2 + xj2 - 2xixj * cosθ), 5)
                    F[i,i] -= (4xi2 + xj2 + xj2 * cos2θ) * d
                    F[i,j] += 4 * (xi2 + xj2) * cosθ * d
                elseif k > 0
                    F[i, j] -= 1 / (sqrt(1 - cosθ) * mypow(x[i], 3))
                end
            end
        end
    end
    return Symmetric(F)
end

```

```julia
julia> all(funS(x, 20) .∈ funS(xi, 20))
true

julia> @btime funS($x, 20);
  2.498 ms (39458 allocations: 630.88 KiB)

julia> @btime funS($xi, 20);
  3.342 ms (46318 allocations: 1.34 MiB)

```

---

<div class="post-metadata">

**Author:** ![cvanaret](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cvanaret/32/11594_2.png) [@cvanaret](https://discourse.julialang.org/u/cvanaret)\
**Post date:** [November 29, 2019, 10:02am UTC](https://discourse.julialang.org/t/interval-arithmetic-computation-time/14633/17 "2019-11-29T10:02:37Z")

</div>

Just a note on using intervals: you can substitute x^2 with x\*x as long as x does not (strictly) contain 0. Otherwise, the [dependency effect](https://en.wikipedia.org/wiki/Interval_arithmetic#Dependency_problem) kicks in and the result will be an overapproximation.