Build faster paths for SAF constructions

[matbesancon]

I was repeatedly solving large LPs that were fast enough and the bottleneck ended up being the SAF constructions. Specifically, for λ a vector of variable indices:

dummy_objective = sum(λ, init=0.0)

was noticeably slow when profiling, and:

nv = length(atoms)
for i in 1:nv
    lhs[i] = 2 * sum(λ[j] * dot(atoms[i], atoms[j]) for j in 1:nv)
    rhs[i] = -dot(atoms[i], b)
end

was also extremely slow (taking more time than a solve. The length of atoms was in the hundreds.

For both, the bulk of the time was spent in the SAF construction, MA.operate, etc

[matbesancon]

Some potential improvements:

1. a MOI manual section on performance gotchas with deceptively slow operations
2. specializing at least common ones? sum and dot for instance
3. ?

[blegat]

We have specialized sum, but you don't hit it because of the keyword argument, that should be fixed in https://github.com/jump-dev/MutableArithmetics.jl/pull/306 You can also use the following, that's what MA will do once that PR is merged:

λ = MOI.VariableIndex.(1:n)
mapreduce(identity, MA.add!!, λ; init = 0.0)

For the second example, use dot, it will use an optimized method in MA:

nv = length(atoms)
for i in 1:nv
    lhs[i] = dot(λ, [2dot(atoms[i], atoms[j]) for j in 1:nv])
    rhs[i] = -dot(atoms[i], b)
end

[matbesancon]

you do need the initial argument though, otherwise it doesn't work for summing variable indices:

julia> sum(l)
ERROR: MethodError: no method matching zero(::Type{MathOptInterface.VariableIndex})
The function `zero` exists, but no method is defined for this combination of argument types.

[matbesancon]

For the second example, use dot, it will use an optimized method in MA:

the only issue is that you need to allocate a whole vector just to construct another SAF which itself will also allocate

[blegat]

Yes, you need the init, that's why we need https://github.com/jump-dev/MutableArithmetics.jl/pull/306

You don't need to allocate anything. dot calls fused_map_reduce which needs getindex to be defined so you cannot use Base.Generator but you can use MOI.Utilities.VectorLazyMap. You can also simply do:

mapreduce(identity, (acc, j) -> add_mul!!(acc, λ[j], 2dot(atoms[i], atoms[j])), 1:nv; init = 0.0)

Or just write the for loop calling add_mul!!