regression observed in JuMP expressions

[kalmarek]

MWE:

using Random
Random.seed!(1)
let dims = [31, 24, 49, 30, 33, 49, 33, 21, 44, 50, 50, 32, 31, 48, 28, 34, 35, 41, 25, 43]
    model = Model()
    Ps = [
        JuMP.@variable(model, [1:d, 1:d])
        for d in dims
    ]
    Ms = [rand(d,d) for d in dims]
    
    @time JuMP.@constraint(model, sum(dot(P, M) for (P,M) in zip(Ps, Ms)) == 1.0 )
    @time JuMP.@constraint(model, sum(dot(Ps[i], Ms[i]) for i in eachindex(dims)) == 1.0 )
    
    model
end;

in JuMP 0.20.1 I get

0.006496 seconds (955 allocations: 5.136 MiB)
0.006837 seconds (873 allocations: 5.133 MiB)

in JuMP 0.21.10

0.034702 seconds (114.34 k allocations: 7.549 MiB)
0.035239 seconds (114.26 k allocations: 7.546 MiB)

which is 5× slowdown (on real examples even larger, closer to 10×.

@blegat

[odow]

Nice catch. It's obviously not going via the mutable route!

[blegat]

I get issues similar to jump-dev/MathOptInterface.jl#1353 where I get allocations when calling a function for no reason I can think of. Copy pasting the same functions outside of the module sometimes decrease the allocation.
See the following gist: https://gist.github.com/blegat/6d7a179e258edcbd31f00443e34d37e4
which gives me:

julia> include("bug.jl");
  0.000660 seconds (3.46 k allocations: 115.484 KiB)
  0.000072 seconds (1.85 k allocations: 77.984 KiB)
  0.000070 seconds (1.85 k allocations: 77.984 KiB)
  0.000087 seconds (2.65 k allocations: 96.734 KiB)
  0.000636 seconds (2.65 k allocations: 96.734 KiB)
  0.000675 seconds (3.46 k allocations: 115.484 KiB)

julia> include("bug.jl");
  0.000685 seconds (3.46 k allocations: 115.578 KiB)
  0.000091 seconds (1.86 k allocations: 78.078 KiB)
  0.000098 seconds (1.86 k allocations: 78.078 KiB)
  0.000114 seconds (2.66 k allocations: 96.828 KiB)
  0.000675 seconds (2.66 k allocations: 96.828 KiB)
  0.000685 seconds (2.66 k allocations: 96.828 KiB)

It seems #116 improves this time-wise (but I still don't understand these allocations). After the PR, I get:

julia> include("bug.jl");
  0.000090 seconds (2.64 k allocations: 96.547 KiB)
  0.000076 seconds (1.84 k allocations: 77.797 KiB)
  0.000080 seconds (1.84 k allocations: 77.797 KiB)
  0.000097 seconds (2.64 k allocations: 96.547 KiB)
  0.000123 seconds (2.64 k allocations: 96.547 KiB)
  0.000110 seconds (2.64 k allocations: 96.547 KiB)

julia> include("bug.jl");
  0.000114 seconds (2.67 k allocations: 96.891 KiB)
  0.000079 seconds (1.86 k allocations: 78.141 KiB)
  0.000082 seconds (1.86 k allocations: 78.141 KiB)
  0.000088 seconds (2.67 k allocations: 96.891 KiB)
  0.000103 seconds (2.67 k allocations: 96.891 KiB)
  0.000107 seconds (2.67 k allocations: 96.891 KiB)

[odow]

Is this another stack/heap heuristic failure?

[blegat]

Is this another stack/heap heuristic failure?

I don't know where it would come from. I use Vararg instead of ... so I'm guessing it's maybe op::Function that should be replaced by op::F ... where {F<:Function}. I don't know of any other possibilities though :/

[blegat]

Is this another stack/heap heuristic failure?

I don't know where it would come from. I use Vararg instead of ... so I'm guessing it's maybe op::Function that should be replaced by op::F ... where {F<:Function} (hence #116 and #120). I don't know

[kalmarek]

@blegat Is this considered to be solved? with

(jl_SlG49P) pkg> st
      Status `/tmp/jl_SlG49P/Project.toml`
  [4076af6c] JuMP v0.22.2

I get

julia> let N = 100
           Random.seed!(1)
           m = JuMP.Model()
           JuMP.@variable(m, P[1:N, 1:N], Symmetric)
           JuMP.@constraint(m, P in JuMP.PSDCone())
           A = sprand(N, N, 0.03)
           B = rand(N, N)
           @time JuMP.@constraint(m, dot(A, P) == 0.0)
           @time JuMP.@constraint(m, dot(B, P) == 0.0)
           m
       end
  0.001497 seconds (80.05 k allocations: 2.947 MiB)
  0.002826 seconds (109.54 k allocations: 3.665 MiB)

also with modified mutable_dot.jl (from gist to MutableArithmetics-0.3.2) I get:

julia> include("mutable_dot.jl");
  0.000052 seconds (817 allocations: 46.953 KiB)
  0.000028 seconds (17 allocations: 28.203 KiB)
  0.000026 seconds (17 allocations: 28.203 KiB)
  0.000042 seconds (817 allocations: 46.953 KiB)
  0.000041 seconds (817 allocations: 46.953 KiB)
  0.000039 seconds (817 allocations: 46.953 KiB)

julia> f(100);
  0.000985 seconds (20.03 k allocations: 1.498 MiB)
  0.000713 seconds (28 allocations: 1.040 MiB)
  0.000784 seconds (28 allocations: 1.040 MiB)
  0.001074 seconds (20.03 k allocations: 1.498 MiB)
  0.001052 seconds (20.03 k allocations: 1.498 MiB)
  0.001047 seconds (20.03 k allocations: 1.498 MiB)

which means that currently

function fused_map_reduce0(
    op::F,
    a, b,
) where {F<:Function,N}
    MA._same_length(a, b)
    T = MA.promote_map_reduce(op, eltype(a), eltype(b))
    accumulator = MA.neutral_element(MA.reduce_op(op), T)
    for I in zip(eachindex(a), eachindex(b))
        JuMP._add_or_set!(accumulator.terms, getindex(a, I[1]), getindex(b, I[2]))
    end
    return accumulator
end

function fused_map_reduce1(
    op::F,
    a, b,
) where {F<:Function,N}
    MA._same_length(a, b)
    T = MA.promote_map_reduce(op, eltype(a), eltype(b))
    accumulator = MA.neutral_element(MA.reduce_op(op), T)
    for I in zip(eachindex(a), eachindex(b))
        accumulator = MA.operate!(
            MA.add_mul,
            accumulator,
            adjoint(getindex(a, I[1])),
            getindex(b, I[2]),
        )
    end
    return accumulator
end

are the fastest

[odow]

Closing in favor of #132