using MutableArithmetic makes DynamicPolynomials slower?

[kalmarek]

Setup: (mt is a multiplication table of an algebra)

using DynamicPolynomials
using MutableArithmetics

@inline function algebra_addmul!(
        res::AbstractVector,
        X::AbstractVector,
        Y::AbstractVector,
        mt::AbstractMatrix{<:Integer})

    @inbounds for j in 1:size(mt,2)
        for i in 1:size(mt,1)
            res[mt[i,j]] = MutableArithmetics.add_mul!(res[mt[i,j]], X[i], Y[j])
#             res[mt[i,j]] += X[i]*X[j] 
        end
    end
    return res
end

function algebra_SOS(Q::AbstractMatrix{T}, mt::AbstractMatrix{<:Integer}) where T
#     result = zeros(T, maximum(mt)) doesn't work
    result = fill(zero(T), maximum(mt))
    for i in 1:size(Q,2)
        @views algebra_addmul!(result, Q[:,i], Q[:,i], mt)
    end
end

and timings (on contrived example):

using Random

let N = 13
    Random.seed!(1)
    @polyvar V[1:N^2]
    mt = rand(1:N^2, N, N)
    Q = V[rand(1:N^2, N, N)]

    algebra_SOS(Q, mt)
    m = rand(N,N)
    algebra_SOS(m, mt)

    @time algebra_SOS(Q, mt)
    @time algebra_SOS(m, mt)
end;

# with res[mt[i,j]] = MutableArithmetics.add_mul!(res[mt[i,j]], X[i], Y[j])
  0.261018 seconds (833.87 k allocations: 47.396 MiB)
  0.000004 seconds (1 allocation: 1.453 KiB)

# with res[mt[i,j]] += X[i]*X[j] 
0.015759 seconds (120.20 k allocations: 15.800 MiB)
  0.000008 seconds (1 allocation: 1.453 KiB)

@blegat, @tweisser

[blegat]

This is weird, MA should be much faster for this with less allocations. Does the issue persist if mt[i, j] is replaced by i + (j - 1) * N ? EDIT it seems it does persist

[blegat]

https://github.com/JuliaAlgebra/DynamicPolynomials.jl/pull/65 solves one part of the performance issue but it's still slower.

[blegat]

The issue is that fill(zero(T), maximum(mt)) uses the same polynomial for every entry so it's slower because it modifies a large polynomial. You should use [zero(T) for i in 1:maximum(mt)] for instance

[kalmarek]

@blegat thanks for investigating; interesting find: indeed v = fill(zero(T), 3) results in identical (===) entries, but modifying one unaliases it from the others (?). When using

function algebra_SOS(Q::AbstractMatrix{T}, mt::AbstractMatrix{<:Integer}) where T
#     result = zeros(T, maximum(mt)) doesn't work
    result = [zero(T) for _ in 1:maximum(mt)]
    for i in 1:size(Q,2)
        @views algebra_addmul!(result, Q[:,i], Q[:,i], mt)
    end
end

I get the following timings:

0.026193 seconds (235.44 k allocations: 30.099 MiB, 17.03% gc time) # with add_mul!
0.015501 seconds (121.20 k allocations: 15.861 MiB) # with res[mt[i,j]] += X[i]*X[j]