performance for expanding large expressions with multilinear functions

[matthias314]

As suggested by @ChrisRackauckas, I want to draw your attention to this Discourse thread about implementing rules for multilinear functions in SymbolicUtils.jl. @greatpet and I came up with a @rule version and one with hand-coded functions. We compared them to a rule-based version in Mathematica and a procedure-based version in Maple. (All code is on Discourse.) Maple and Mathematica seem to be comparable, with Maple being somewhat faster. (Disclaimer: I don't have access to Mathematica myself.)

It turned out that the rule-based version in SymbolicUtils.jl is much slower and not usable on bigger examples. The hand-coded version for a multilinear function f is comparable to Maple for non-trivial examples like f(a,a,a,a,a) with a=u+2v+3w+4x+5y+6z (which leads to close to 8000 terms). However, on larger examples performance drops. With 7 arguments equal to a (close to 280,000 terms), it's 20 times slower than Maple. Maybe there is room to improve the performance of your package -- or maybe we were simply not doing things the right way. (For example, would performance improve if we used Symbolics.jl instead?)

EDIT: Here are some timings (in seconds) of the Julia 1.6.1/SymbolicUtils 0.13.1 rule and code versus Maple 2019. (I've used a fairly old machine on which Maple is installed.) f and a are as above. n is the number of a arguments to f. Remember that the Maple time() function is not very precise.

n	Maple	Julia code	Julia rule
4	0.032	0.014	2.0
5	0.19	0.21	20.6
6	1.15	4.97
7	8.8	220

For the n = 7 computation Maple uses 137 MB and Julia maybe 480 MB.

[zengmao]

For completeness, here's the rule-based Mathematica code, which was not fully written out in the Discourse thread,

rule = {
    f[a___, HoldPattern[Plus[b__]], c___] :>
        Module[{i},
            Total[Table[f[a, List[b][[i]], c], {i, Length[List[b]]}]]
        ],
    f[a___, coeff_Integer * b_, c___] :> coeff * f[a,b,c]
};

Timing[
f[u+2v+3w+4x+5y+6z, u+2v+3w+4x+5y+6z, u+2v+3w+4x+5y+6z, u+2v+3w+4x+5y+6z, u+2v+3w+4x+5y+6z] //. rule;
]

Below is a comparison of timing between Mathematica, the "imperative" Julia code, and the rule-based Julia code on my machine. (My Maple license expired recently.) You can see that the hand-coded Julia version is the fastest here, but when you write rule-based code, SymbolicUtils.jl is currently much slower than Mathematica.

n	Mathematica	Julia code	Julia rule
4	0.019	0.0075	0.90
5	0.15	0.12	11

In case it's handy, the rule-based Julia version is copied here, with timing for n=5 terms.

using SymbolicUtils, SymbolicUtils.Rewriters

using BenchmarkTools

const f = SymbolicUtils.Sym{SymbolicUtils.FnType{Tuple, Number}}(:f)

r1 = @rule(f(~~a, +(~~b), ~~c) => sum([f(~~a..., x, ~~c...) for x in ~~b]))
r2 = @rule(f(~~a, ~c::(x -> x isa Int) * ~x, ~~b) => ~c * f(~~a..., ~x, ~~b...))

expand_multilinear = Fixpoint(Postwalk(Chain([r1, r2])))

@syms u v w x y z

@btime expand_multilinear(f(u+2v+3w+4x+5y+6z, u+2v+3w+4x+5y+6z, u+2v+3w+4x+5y+6z, u+2v+3w+4x+5y+6z, u+2v+3w+4x+5y+6z));