Extra allocations from multiplication (redefining functions in the REPL alleviates allocations)

[OlivierHnt]

This is a follow-up of this discussion in discourse.

Multiplication of intervals in functions defined inside a module leads to 6 extra allocations per *. Interestingly, redefining the function in the REPL fixes the issue.

Here is a MWE:

module MyModule
    foo(x) = x*x
    export foo
end

using .MyModule, IntervalArithmetic, BenchmarkTools

x = @interval(π);

@btime foo(Ref($x)[]);
# 137.597 ns (6 allocations: 96 bytes)

MyModule.foo(x) = x*x

@btime foo(Ref($x)[]);
# 59.040 ns (0 allocations: 0 bytes)

[OlivierHnt]

To add to this issue: the division operator / does not have this problem and * has the extra allocations only when multiplying two intervals.

See the MWE:

module MyModule
    mymul(x, y) = x*y
    mydiv(x, y) = x/y
    export mymul, mydiv
end

using .MyModule, IntervalArithmetic, BenchmarkTools

x = @interval(π);

y = @interval(π);
@btime mymul($x, $y);
@btime mydiv($x, $y);

y = π;
@btime mymul($x, $y);
@btime mydiv($x, $y);

I am not familiar with the internal choices of the package, but could * be reworked in a similar style to /?

[Kolaru]

It is because somehow inlining of foo doesn't happen. If you explicitly @inline it (as it is done internally by IntervalArithmetic for *) the extra allocations are gone

julia> module CoolModule
           foo(x) = x*x
           @inline infoo(x) = x*x
           export foo, infoo
       end
Main.CoolModule

julia> using IntervalArithmetic

julia> using BenchmarkTools

julia> using .CoolModule

julia> x = @interval(π)
[3.14159, 3.1416]

julia> @btime foo(Ref($x)[])
  43.116 ns (6 allocations: 96 bytes)
[9.8696, 9.86961]

julia> @btime infoo(Ref($x)[])
  10.372 ns (0 allocations: 0 bytes)
[9.8696, 9.86961]

We have to specifically inline * and not / because the former has some annoying corner cases, leading to a more involved code the compiler may be reluctant to inline.

That being sayed the issue remains: it is not good to loose the inlining benefit in user defined modules.

[dpsanders]

Thanks for the detective work!

We should probably just use @inline on all these functions.

[OlivierHnt]

As I mentioned on discourse, unfortunately the macro @inline does not always prevent the extra allocations.

Here is yet an other behaviour:

module MyModule
    using IntervalArithmetic
    *(Interval(1), Interval(1)) # <- prevents extra allocations

    function foo(v)
        w = zero.(v)
        for i ∈ eachindex(w)
            _foo!(w, v, i)
        end
        return w
    end

    function _foo!(w, v, i)
        for j ∈ 1:i
            w[i] += v[i]*v[j]
        end
    end
   
    export foo
end

using .MyModule, IntervalArithmetic, BenchmarkTools

v = Interval.(rand(1_000));
@btime foo($v);

It is a convenient work around, however it does not seem to scale well with the project size. I tried it on a module which include multiple files and it does not work; for large projects, re-defining the functions remain the only way to prevent allocations as far as I know.

[OlivierHnt]

The following operation is very slow and may be related to this issue:

using IntervalArithmetic, BenchmarkTools

A = Interval.(rand(100, 100))

@btime $A * $A; # 184.122 ms (6000020 allocations: 91.71 MiB) 

It is becoming less and less easy to hack around this issue 😞. (See edit below)

I tried to look more into it (e.g. with some profiling)... Not much luck.

@Kolaru have you been able to investigate further?
Also I wonder if this problem also occurs with TaylorModels.jl.

EDIT: actually there seems to be a convenient hack, it seems that one can just redefine the problematic * in the REPL after compilation, that is to say:

using IntervalArithmetic

#-- prevents allocations
Interval(1.0) * Interval(1.0)

@inline function Base.:*(a::Interval{T}, b::Interval{T}) where T<:Real
    (isempty(a) || isempty(b)) && return emptyinterval(T)

    (iszero(a) || iszero(b)) && return zero(Interval{T})

    (isfinite(a) && isfinite(b)) && return IntervalArithmetic.mult(*, a, b)

    return IntervalArithmetic.mult(IntervalArithmetic.checked_mult, a, b)
end
#---

# the real code starts here...

[dpsanders]

You can also use setrounding(Interval, :fast) to use a faster method of directed rounding.

If you really need fast multiplication of interval matrices then you will need to implement Rump's method using rounding-mode changes and floating-point matrix multiplication (can't find the reference right now).

[OlivierHnt]

Thank you for your feedback. Are you referring to this paper Fast interval matrix multiplication?

[lucaferranti]

You can find the documentation of different rounding modes here.

[OlivierHnt]

Thanks for the link.

I went through the five rounding modes:

using IntervalArithmetic, BenchmarkTools

A = Interval.(rand(100, 100));

setrounding(Interval, :fast);
@btime $A*$A; # 69.620 ms (8 allocations: 156.66 KiB)

setrounding(Interval, :accurate);
@btime $A*$A; # 13.685 ms (8 allocations: 156.66 KiB)

setrounding(Interval, :tight);
@btime $A*$A; # 172.094 ms (6000020 allocations: 91.71 MiB)

setrounding(Interval, :slow);
@btime $A*$A; # 226.559 ms (8 allocations: 156.66 KiB)

setrounding(Interval, :none);
@btime $A*$A; # 8.206 ms (8 allocations: 156.66 KiB)

#-- trick to remove extra allocations observed when rounding set to :tight
@inline function Base.:*(a::Interval{T}, b::Interval{T}) where T<:Real
    (isempty(a) || isempty(b)) && return emptyinterval(T)

    (iszero(a) || iszero(b)) && return zero(Interval{T})

    (isfinite(a) && isfinite(b)) && return IntervalArithmetic.mult(*, a, b)

    return IntervalArithmetic.mult(IntervalArithmetic.checked_mult, a, b)
end
#---

setrounding(Interval, :tight);
@btime $A*$A; # 84.958 ms (8 allocations: 156.66 KiB)

The extra allocations are only observed when the rounding is set to :tight; perhaps you already knew this though :).

[dpsanders]

@OlivierHnt: Have you tried inlining the functions inside the package?

[OlivierHnt]

Yes I did try it on a large project but it did not work.
Somehow, sometimes (e.g. for setrounding(Interval, :tight)), the inlining of * is lost.

Additionally, in the MWE @Kolaru provided with @inline there was no need to call using IntervalArithmetic inside the module. I gave a counter example on Discourse: on a slightly bigger module, it becomes necessary (cf. the following code).

Code

using IntervalArithmetic, BenchmarkTools

v = Interval.(rand(1_000));

module M1    
    function M1_foo(v)
        w = zero.(v)
        for i ∈ eachindex(w)
            _foo!(w, v, i)
        end
        return w
    end

    @inline function _foo!(w, v, i)
        for j ∈ 1:i
            w[i] += v[i]*v[j]
        end
    end
   
    export M1_foo
end

using .M1
@btime M1_foo($v);

module M2
    using IntervalArithmetic
    
    function M2_foo(v)
        w = zero.(v)
        for i ∈ eachindex(w)
            _foo!(w, v, i)
        end
        return w
    end

    @inline function _foo!(w, v, i)
        for j ∈ 1:i
            w[i] += v[i]*v[j]
        end
    end
   
    export M2_foo
end

using .M2
@btime M2_foo($v);

Also, these sort of things happen:

using IntervalArithmetic, BenchmarkTools

A = Interval.(rand(10, 10));

Interval(1)*Interval(1) # we perform `*` once before any computation

@btime $A*$A; # 85.819 μs (5 allocations: 1.98 KiB)

compared to

using IntervalArithmetic, BenchmarkTools

A = Interval.(rand(10, 10));

A*A # some computation...

@btime $A*$A; # 191.300 μs (6017 allocations: 95.92 KiB)

Interval(1)*Interval(1) # this will no longer prevent the extra allocations

@btime $A*$A; # 191.681 μs (6017 allocations: 95.92 KiB)

If we give Julia the opportunity to compile just *, things seem fine. However, if we perform an other operation before, then Julia seems to get mixed up...

[dpsanders]

I meant: Did you try modifying the code in IntervalArithmetic.jl to explicitly inline all the relevant functions?

[OlivierHnt]

I see. No I did not try it; I'll give it a shot perhaps next week then.

[OlivierHnt]

Sorry for the late feedback. I used @inline in front of the following functions:

• IntervalArithmetic.jl/src/intervals/rounding.jl

  Lines 134 to 135 in a3dc27f

  	@eval $op(::IntervalRounding{:tight},
  	a::$T, b::$T, $mode1) = $ff(a, b)

• IntervalArithmetic.jl/src/intervals/rounding.jl

  Lines 175 to 180 in a3dc27f

  	@eval function $f(::IntervalRounding{:tight},
  	a::T, b::T, $mode1) where T<:AbstractFloat
  	setrounding(T, $mode2) do
  	$f(a, b)
  	end
  	end

• IntervalArithmetic.jl/src/intervals/rounding.jl

  Line 278 in a3dc27f

  @eval $f(a::T, b::T, r::RoundingMode) where {T<:AbstractFloat} = $f($roundtype, a, b, r)

No luck.

[dpsanders]

@OlivierHnt: What about mult for example? And any other function involved in the multiplication.