ScalarFunctionIterator has quadratic complexity

[mohamed82008]

The current ScalarFunctionIterator has a best case quadratic and worst case cubic (if the constraints are dense) complexity. If we do a Radix sort by the output_index then it becomes O(nk) where k is the average number of terms per constraint. For sparse constraints, this is linear complexity. This together with the use of StructsOfArrays.jl can enable returning a view from scalar_terms_at_index which reduces allocations. The catch here is if the user/solver expects the terms to be in a certain order and sorting will just confuse everyone. I haven't seen this assumption made anywhere so I am asking if you know otherwise.

https://github.com/JuliaOpt/MathOptInterface.jl/blob/247c4baf673eb8b1c29177e7f1652a0a78e1d570/src/Utilities/functions.jl#L103

[mohamed82008]

Perhaps the StructsOfArrays point is a discussion for another thread since it involves restructuring not just an implementation detail.

[odow]

We should hold off until we can benchmark the performance. Where is this used in the rest of the code?

[mohamed82008]

Where is this used in the rest of the code?

I am still reading so I am not sure but this seems like the kind of thing I would assume is efficient since it is baked in MOI. Maybe others know where it is used.

[mohamed82008]

Shouldn't be hard to benchmark this specific part of the code using random sparse matrices of coefficients converted to Vector{VectorAffineTerm} which are then converted to VAFs, then iterated over. I will try to do it and see what happens.

[blegat]

getindex used to be called directly on the AbstractVectorFunction.
Now that it eachscalar is called first, we could add caching in the ScalarFunctionIterator to avoid having to redo findall for every getindex.
We could for instance have a Vector{Vector{Int}} stored in cache storing the list of indices of terms for a specific output index. This can be constructed from a single pass through the terms in eachscalar.

[mohamed82008]

This is a nice middle way solution that does not involve the hidden sorting of the terms. Sorting will probably use less memory but may be an unexpected and unwelcome side effect.

[blegat]

Yes, we need to consider the function as read-only :-P

[mohamed82008]

If memory is not of primary concern, we can also cache a sorted version of the vector of terms and a vector of end markers for the rows' terms. This will enable using the StructsOfArrays trick, which will save allocations by returning a cache-friendly view into the sorted vector of terms.

Either way, we can return a view, but the sorted one is more cache friendly since it is all contiguous in memory.

[blegat]

The ScalarAffineFunction terms need to be a Vector, not any AbstractVector so it cannot be a view, can it ?

[mohamed82008]

True, but that's something that will also have to change :) Parameterizing by the vector type is the first step towards native GPU solver support. It's also good for cases like this to avoid allocations using views.

[mohamed82008]

That's a discussion for another thread though.

[blegat]

This will completely copy all the terms of the function in the sorting so if the user only wants to obtain the scalar function at one index, it might be overkill to create the list of sorted terms.
Putting the indices of the terms of each output index in different buckets is the first step of the radix sort so it makes sense to cache that and have fast creation of scalar function when needed.
If we parametrize by the vector type the view idea makes sense. IIRC, it has already been discussed and we decided against the idea to keep a simpler interface. We might reconsider but it might take time to change this. If you want this issue to be fixed in a reasonable amount of time, it is probably best not to rely on this :)

[mohamed82008]

Coming up in a PR.