`MOI.Utilities.operate` StackOverflowError

[hdavid16]

The following returns a StackOverflowError due to deep recursion:

x = [i + i*MOI.VariableIndex(i) for i in 1:1000]
MOI.Utilities.operate(+, Int, x...)

operate is used in SolverAPI.jl to convert from SNF to SAF or SQF when possible.

The proposed fix is to rewrite the following

MathOptInterface.jl/src/Utilities/operate.jl

Lines 323 to 327 in 53ec269

	### 1c: operate(+, T, args...)
	function operate(::typeof(+), ::Type{T}, f, g, h, args...) where {T<:Number}
	return operate!(+, T, operate(+, T, f, g), h, args...)
	end

as (credit to @bachdavi for this solution):

function operate(::typeof(+), ::Type{T}, args...) where {T<:Number}
    operate_plus(accum, x) = MOI.Utilities.operate!(+, T, accum, x)
    return reduce(operate_plus, args)
end

Some tests shown below (Note: operate2 is the new form of operate above):

# 100 ScalarAffineFunction args
julia> x = [i + i*MOI.VariableIndex(i) for i in 1:100];

julia> @time MOIU.operate(+, Int, x...);
  0.000200 seconds (199 allocations: 46.078 KiB)

julia> @time MOIU.operate2(+, Int, x...);
  0.000030 seconds (7 allocations: 6.797 KiB)
#--------------------------------------------------
# 1000 ScalarAffineFunction{Int} args  
julia> x = [i + i*MOI.VariableIndex(i) for i in 1:1000];

julia> @time MOIU.operate(+, Int, x...);
ERROR: StackOverflowError:

julia> @time MOIU.operate2(+, Int, x...);
  0.000093 seconds (9 allocations: 74.734 KiB) 
#--------------------------------------------------
# 1000 ScalarAffineFunction{Float64} args
julia> x = [i + i*MOI.VariableIndex(i) for i in 1.0:1000.0];

julia> @time MOIU.operate2(+, Float64, x...);
  0.000116 seconds (9 allocations: 74.734 KiB)
#--------------------------------------------------  
# 1000 ScalarQuadraticFunction args
julia> x = [i + i*MOI.VariableIndex(i)*MOI.VariableIndex(i) for i in 1:1000];

julia> @time MOIU.operate(+, Int, x...);
ERROR: StackOverflowError:

julia> @time MOIU.operate2(+, Int, x...);
  0.000105 seconds (5 allocations: 127.562 KiB)
#--------------------------------------------------  
# 1000 mixed aff and quad args
julia> x = [i + i*MOI.VariableIndex(i)*MOI.VariableIndex(i) for i in 1:500];

julia> y = [i + i*MOI.VariableIndex(i) for i in 1:500];

julia> z = vcat(x,y);

julia> @time MOIU.operate(+, Int, z...);
ERROR: StackOverflowError:

julia> @time MOIU.operate2(+, Int, z...);
  0.000126 seconds (6 allocations: 102.750 KiB)

Moving to higher number of affine function args:

Happy to submit a PR to MOI if desired.

[odow]

I don't know if we every considered splatting with such a high number of args. You can do instead operate(sum, x). I thought this was MA.operate

[odow]

We can fix this particular case, but it hints at a larger problem in your code.

[odow]

Fixed in #2285

For well-typed arrays, Base.sum is faster:

julia> @time MOI.Utilities.operate(+, Int, x...);
  0.000041 seconds (8 allocations: 51.281 KiB)

julia> @time sum(x);
  0.000074 seconds (7 allocations: 43.328 KiB)

But it doesn't work for abstractly-typed arrays

julia> x = [i + i*MOI.VariableIndex(i)*MOI.VariableIndex(i) for i in 1:500];

julia> y = [i + i*MOI.VariableIndex(i) for i in 1:500];

julia> z = vcat(x,y);

julia> @time sum(z);
ERROR: MethodError: no method matching zero(::Type{MathOptInterface.AbstractScalarFunction})

[hdavid16]

Thanks. Yes, Base.sum works well in some cases. Since we are taking the args form a ScalarNonlinearFunction and attempting to convert it into a ScalarAffineFunction or ScalarQuadraticFunction, we have a Vector{Any}, which works fine with Base.sum, unless the args are all MOI.VariableIndex.

[blegat]

I don't agree with #2285, it's going to be type unstable to do a for loop in a tuple that contains elements of different types.
I think it makes sense for this to be a StackOverflow, you shouldn't use splatting when you have that many arguments.
If you have a mix of ScalarAffineFunction and ScalarQuadraticFunction, you can either convert everything to ScalarQuadraticFunction or have one vector of ScalarAffineFunction and one of ScalarQuadraticFunction and sum one and then the other one. Or you can write your own sum. Doing the for loop of #2285 on a Vector{Any} (or is it a tuple ?) is probably not going to be fast anyway I guess.

[odow]

I don't agree with #2285, it's going to be type unstable to do a for loop in a tuple that contains elements of different types.
I think it makes sense for this to be a StackOverflow, you shouldn't use splatting when you have that many arguments.

I pretty strongly disagree. We've spent too much effort focusing on performance. This was clearly a bug because a user encountered it in the wild. There's probably a better way to write it that doesn't involve a 1,000 argument splat, but we shouldn't stackoverflow.

[odow]

Do you have an example where this is a performance problem?

[blegat]

I pretty strongly disagree. We've spent too much effort focusing on performance. This was clearly a bug because a user encountered it in the wild. There's probably a better way to write it that doesn't involve a 1,000 argument splat, but we shouldn't stackoverflow.

We use recursion on the tuple of arguments pretty much everywhere in the code, are we going to rewrite everything with for loops ?
I think the recursion is fine, if we want to also make it work well with a long number of arguments, we can use afoldl (like + does for instance) but I don't agree with using the for loops all the time.
It's pretty easy to build an example, if you use MA.operate!(+, ::SAF, ::VI, ::SAF, ::VI, ::SAF), you get almost no allocation with the current approach and #2285 will be type unstable hence will allocate.

[odow]

It's pretty easy to build an example, if you use MA.operate!(+, ::SAF, ::VI, ::SAF, ::VI, ::SAF), you get almost no allocation with the current approach and #2285 will be type unstable hence will allocate.

This is not true. My way is faster and allocates less, although it is type unstable

(moi) pkg> st
Status `/private/tmp/moi/Project.toml`
  [b8f27783] MathOptInterface v1.20.0

julia> import MathOptInterface as MOI

julia> using BenchmarkTools

julia> function my_operate(::typeof(+), ::Type{T}, f, g, h, args...) where {T<:Number}
           ret = MOI.Utilities.operate(+, T, f, g)
           ret = MOI.Utilities.operate!(+, T, ret, h)
           for a in args
               ret = MOI.Utilities.operate!(+, T, ret, a)
           end
           return ret
       end
my_operate (generic function with 1 method)

julia> f = (
           MOI.VariableIndex(1),
           1 * MOI.VariableIndex(1) + 2,
           3 * MOI.VariableIndex(1) * MOI.VariableIndex(1) + 4,
       )
(MOI.VariableIndex(1), (2) + (1) MOI.VariableIndex(1), (4) + 3.0 MOI.VariableIndex(1)²)

julia> x = f[[2, 1, 2, 1, 2, 3]]
((2) + (1) MOI.VariableIndex(1), MOI.VariableIndex(1), (2) + (1) MOI.VariableIndex(1), MOI.VariableIndex(1), (2) + (1) MOI.VariableIndex(1), (4) + 3.0 MOI.VariableIndex(1)²)

julia> foo(x) = MOI.Utilities.operate(+, Int, x...)
foo (generic function with 1 method)

julia> bar(x) = my_operate(+, Int, x...)
bar (generic function with 1 method)

julia> @benchmark foo($x)
BenchmarkTools.Trial: 10000 samples with 151 evaluations.
 Range (min … max):  679.543 ns …  21.893 μs  ┊ GC (min … max): 0.00% … 95.60%
 Time  (median):     699.675 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   757.773 ns ± 810.786 ns  ┊ GC (mean ± σ):  4.47% ±  4.05%

  ▃█▆▃▁      ▁           ▁                                      ▁
  ██████▆▇▇▇▇█▆█▇█▅▅▅▇▇▇▇█▇▆▇▅▄▅▅▅▃▄▄▃▄▅▅▅▄▅▅▄▄▅▅▁▅▆▁▅▅▅▅▄▅▄▄▅▆ █
  680 ns        Histogram: log(frequency) by time       1.39 μs <

 Memory estimate: 672 bytes, allocs estimate: 11.

julia> @benchmark bar($x)
BenchmarkTools.Trial: 10000 samples with 157 evaluations.
 Range (min … max):  657.471 ns …  21.416 μs  ┊ GC (min … max): 0.00% … 95.95%
 Time  (median):     680.608 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   729.359 ns ± 845.134 ns  ┊ GC (mean ± σ):  4.87% ±  4.07%

   ▂▆██▆▅▃▁                                                     ▂
  ▇████████▇▇▇▄▅▅▄▃▅▄▄▄▆▆▅▄▄▆▅▄▄▆▅▄▄▄▅▅▃▄▆▇▆▇▆▆▇▆▆▄▅▆▄▆▅▄▃▄▂▂▄▃ █
  657 ns        Histogram: log(frequency) by time       1.01 μs <

 Memory estimate: 608 bytes, allocs estimate: 9.

julia> @code_warntype foo(x)
MethodInstance for foo(::Tuple{MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.VariableIndex, MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.VariableIndex, MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.ScalarQuadraticFunction{Int64}})
  from foo(x) @ Main REPL[40]:1
Arguments
  #self#::Core.Const(foo)
  x::Tuple{MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.VariableIndex, MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.VariableIndex, MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.ScalarQuadraticFunction{Int64}}
Body::MathOptInterface.ScalarQuadraticFunction{Int64}
1 ─ %1 = MathOptInterface.Utilities::Core.Const(MathOptInterface.Utilities)
│   %2 = Base.getproperty(%1, :operate)::Core.Const(MathOptInterface.Utilities.operate)
│   %3 = Core.tuple(Main.:+, Main.Int)::Core.Const((+, Int64))
│   %4 = Core._apply_iterate(Base.iterate, %2, %3, x)::MathOptInterface.ScalarQuadraticFunction{Int64}
└──      return %4


julia> @code_warntype bar(x)
MethodInstance for bar(::Tuple{MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.VariableIndex, MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.VariableIndex, MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.ScalarQuadraticFunction{Int64}})
  from bar(x) @ Main REPL[41]:1
Arguments
  #self#::Core.Const(bar)
  x::Tuple{MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.VariableIndex, MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.VariableIndex, MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.ScalarQuadraticFunction{Int64}}
Body::Union{MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.ScalarQuadraticFunction{Int64}}
1 ─ %1 = Main.my_operate::Core.Const(my_operate)
│   %2 = Core.tuple(Main.:+, Main.Int)::Core.Const((+, Int64))
│   %3 = Core._apply_iterate(Base.iterate, %1, %2, x)::Union{MathOptInterface.ScalarAffineFunction{Int64}, MathOptInterface.ScalarQuadraticFunction{Int64}}
└──      return %3

To prove that foo is actually calling the right method:

julia> @which MOI.Utilities.operate(+, Int, x...)
operate(::typeof(+), ::Type{T}, f, g, h, args...) where T<:Number
     @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/LQvlf/src/Utilities/operate.jl:325

MathOptInterface.jl/src/Utilities/operate.jl

Lines 325 to 327 in 82e5e58

	function operate(::typeof(+), ::Type{T}, f, g, h, args...) where {T<:Number}
	return operate!(+, T, operate(+, T, f, g), h, args...)
	end

This is a bug fix, and I don't think we should revert the PR.