Slowdown with build_function output in recent versions

[MarcBerliner]

I recently updated to the latest ModelingToolkit version and I've noticed a pretty sizable slowdown in the in-place result from build_function. One 37-line function that I made is 5x slower and a 148-line function is 20x slower compared to using v4.5.3 (I can DM them to the devs). Much longer functions have a more modest slowdown, but it is still quite noticeable. I believe it's due to the change in syntax in recent versions.

(Apologies if this should be for SymbolicUtils.)

[ChrisRackauckas]

Can you share it with me @YingboMa and @shashi ? I was told that all of the specializations were the same in the newest version, but apparently one must have been missed.

[MarcBerliner]

I didn't realize you can't DM on github so I will email the files to you.

[shashi]

Hi @MarcBerliner could you try out this branch #801 and see if it fixes the regression?

Thanks!

[MarcBerliner]

Hi @shashi, I tried that branch but it didn't seem to fix the issue. The f function is about the same, but the g function is even slower because var"##out#364"[CartesianIndex(i, j)] = 0 is littered at every non-zero entry in the matrix.

[ChrisRackauckas]

@MarcBerliner try setting skipzeros=true in build_function. @shashi I think you accidentally changed the default there.

[MarcBerliner]

@ChrisRackauckas setting skipzeros=true fixes that problem, but it's still slower than before, 50 microseconds compared to 40 microseconds with 5.6.3. It's possible the issue is updating CartesianIndex(i, j) instead of nzval of the sparse matrix?

[shashi]

In the code you sent me it does use nzval though!

[shashi]

The f function is about the same

Do you mean about the same between v4 and v5? Or that it is about as slow on v5 as before?

[MarcBerliner]

The #801 branch uses CartesianIndex(i, j), 5.6.3 uses nzval. I mean that using branch #801 is as slow as v5.

[YingboMa]

Are you sure you have deved MTK and using the right branch? That branch doesn't use CartesianIndex at all.

[MarcBerliner]

I was using a clone of the ddcba6c repo. I just tried 7d5839d and build_function now fails for my large sparse matrix with the error message:

ERROR: MethodError: no method matching +(::AbstractAlgebra.Generic.MPoly{BigInt})
Closest candidates are:
  +(::Union{AbstractAlgebra.RingElem, AbstractFloat, Integer, Rational}, ::AbstractAlgebra.Generic.LaurentPolyWrap) at /Users/marc/.julia/packages/AbstractAlgebra/aF5Iw/src/generic/LaurentPoly.jl:133
  +(::T, ::AbstractAlgebra.ResElem{T}) where T<:AbstractAlgebra.RingElem at /Users/marc/.julia/packages/AbstractAlgebra/aF5Iw/src/generic/Residue.jl:308
  +(::T, ::AbstractAlgebra.ResFieldElem{T}) where T<:AbstractAlgebra.RingElem at /Users/marc/.julia/packages/AbstractAlgebra/aF5Iw/src/generic/ResidueField.jl:313
  ...
Stacktrace:
 [1] substitute(::Term{Real}, ::Dict{Sym,AbstractAlgebra.Generic.MPoly}; fold::Bool) at /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/api.jl:51

It still seems to work for a toy problem though. Here's an example showing it using CartesianIndex on 7d5839d:

julia> using ModelingToolkit, SparseArrays

(@v1.5) pkg> status ModelingToolkit
Status `~/.julia/environments/v1.5/Project.toml`
  [961ee093] ModelingToolkit v5.6.3 `https://github.com/SciML/ModelingToolkit.jl.git#7d5839d347`

julia> @variables x
(x,)

julia> y = sparse(1:3,1:3,x)
3×3 SparseMatrixCSC{Num,Int64} with 3 stored entries:
  [1, 1]  =  x
  [2, 2]  =  x
  [3, 3]  =  x

julia> build_function(y,x, skipzeros=true)[2]
:(function (var"##out#461", x)
      #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:252 =#
      #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:253 =#
      begin
          #= /Users/marc/.julia/packages/ModelingToolkit/Mpmku/src/build_function.jl:302 =#
          #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:299 =# @inbounds begin
                  #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:295 =#
                  var"##out#461"[CartesianIndex(1, 1)] = x
                  var"##out#461"[CartesianIndex(2, 2)] = x
                  var"##out#461"[CartesianIndex(3, 3)] = x
                  #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:297 =#
                  nothing
              end
      end
  end)

Compared to using v5.6.3:

julia> using ModelingToolkit, SparseArrays

(@v1.5) pkg> status ModelingToolkit
Status `~/.julia/environments/v1.5/Project.toml`
  [961ee093] ModelingToolkit v5.6.3

julia> @variables x
(x,)

julia> y = sparse(1:3,1:3,x)
3×3 SparseMatrixCSC{Num,Int64} with 3 stored entries:
  [1, 1]  =  x
  [2, 2]  =  x
  [3, 3]  =  x

julia> build_function(y,x, skipzeros=true)[2]
:(function (var"##out#253", x)
      #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:252 =#
      #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:253 =#
      begin
          #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:295 =#
          (var"##out#253").nzval[1] = x
          (var"##out#253").nzval[2] = x
          (var"##out#253").nzval[3] = x
          #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:297 =#
          nothing
      end
  end)

If I'm doing something incorrectly please let me know.

[ChrisRackauckas]

Did this get cleared up?

[isaacsas]

@shashi’s updates last week made a big difference, but I was still seeing like a 50-60% decrease in the generated ODE rhs performance. (With BCR.)