I’m looking for a Julia package that can rewrite symbolic expressions involving functions that don’t have a fixed number of arguments. For example, it should be easy to declare that some function f is multilinear (= linear in each argument if all other arguments are kept fixed) without having to create a separate transformation rule for each possible case (“linear in the k-th argument out of n”).
The package shouldn’t call Python under the hood. I once tried out SymPy as well as pure Python and found them way too slow. Any suggestions?
1 Like
SymbolicUtils.jl, which underlies Symbolics.jl, has a rewrite engine, but I cannot find support for sequence patterns in the documentation (here the sequence would be the list of arguments). This kind of task is easy using sequence patterns in Wolfram Mathematica, but I’d love to learn about a Julia solution if it exists! (I’m trying to move some of my computing tasks from Mathematica to Julia.) Here’s the Mathematica solution, just to make it clear what I meant by sequence patterns etc.
Input:
f[arguments1___, coeff_?NumericQ * variable_, arguments2___] :=
coeff * f[arguments1, variable, arguments2];
f[x, 2*y]
Output:
2 f[x, y]
Input:
f[x, 2*y, 3*z, w]
Output:
Out[3]= 6 f[x, y, z, w]
Just use Symbolics.jl with arrays.
1 Like
Can you give an example? For instance, how would the rule that @greatpet has given in Mathematica above look like in Symbolics.jl? I would be interested in that because it’s exactly of the kind I’m looking for.
1 Like
I’ve cooked up a solution using Symbolics.jl. The function is here, and some example usage follows,
function expand_multilinear(expr)
# check if expr is a single term before proceeding
if !isa(expr, Num) || !isa(expr.val, SymbolicUtils.Term)
return expr
end
function_name = expr.val.f
function_arguments = deepcopy(expr.val.arguments)
overall_coeff = 1
for i in 1:length(function_arguments)
# check if this function argument is a product of things
if isa(function_arguments[i], SymbolicUtils.Mul)
current_coeff = function_arguments[i].coeff
# absorb nontrivial coefficient into an overall factor
if current_coeff != 1
overall_coeff *= current_coeff
function_arguments[i] = function_arguments[i] / current_coeff
end
end
end
overall_coeff * function_name(function_arguments...)
end
Example usage:
julia> using Symbolics
julia> @variables x,y,z,g(..);
julia> expand_multilinear(g(x, 2*y, 3*z))
6g(x, y, z)
julia> expand_multilinear(g(1/x, 2*y*z))
2g(x^-1, y*z)
The main limitation is that this only works for a single term. It gives up if you run e.g. expand_multilinear(g(2*x,y) + g(3*z, y)) . If you have a general expression, presumably you can look at all sub-expressions (i.e. visiting tree nodes of the interal representation) and apply the above function, but I know far too little about Symbolics.jl to figure out how to do this efficiently.
1 Like
Thanks a lot! Although the code doesn’t quite reach the elegance of Mathematica’s rule-based approach … 
It seems that such patters exist under the name of segment patterns, see this part of the documentation. Example:
using SymbolicUtils
@syms x y z
r = @rule(+(~x,~~y) => +(~~y...))
r(x+y+z) # gives y+z
This could be the basis of a very clean solution without writing as much code as in the function expand_multilinear above. However, while arithmetic operators like + accept a variable number of arguments, I don’t know if it’s possible to define a new symbolic function with this property. Does anybody else know?
UPDATE: The following seems to work:
using SymbolicUtils
@syms w x y z
f = SymbolicUtils.Sym{(SymbolicUtils.FnType){Tuple, Number}}(:f)
r = @rule(f(~~a, ~x + ~y, ~~b) => f(~~a..., ~x, ~~b...) + f(~~a..., ~y, ~~b...))
r( f(x+y) ) # f(x) + f(y)
r( f(x+y, z) ) # f(x, z) + f(y, z)
r( f(w, x+y, z) ) # f(w, x, z) + f(w, y, z)
The crucial point is that by defining f without the help of @syms, one can specify a general Tuple type instead of one with a specific number of parameters (as in Tuple{Number,Number}). (I’ve got this idea from the @variables f(..) macro of Symbolics.jl.)
3 Likes
Nicely done, but on further testing, your code returns nothing when you call it with 3 variables added together in an argument slot, like r( f(w, x+y+z, z) ). I haven’t been able to fix it, but here’s an attempt:
r1 = @rule(f(~~a, ~x + ~~y, ~~b) => f(~~a..., ~x, ~~b...) + f(~~a..., +(~~y...), ~~b...))
r1(f(w, x+y+z, z)) # returns f(w, x, z) + f(w, y + z, z)
Do you known a good way to apply the rule recursively to finish up the job?
An unrelated question / feature request for Symbolics.jl: can associative-commutative pattern matching, i.e. @acrule, be applied to arbitrary user-defined associative / commutative functions, rather than Base.+ and Base.*?
Yes, it was not a full solution. Here is what I have at the moment:
using SymbolicUtils, SymbolicUtils.Rewriters
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...))
r = Fixpoint(Postwalk(Chain([r1, r2])))
I haven’t tested it fully, but it seems to work:
julia> @syms w x y z;
julia> r( f(2*w, x-y+z, f(w-3*z,x)) )
2f(w, x, f(w, x)) + 6f(w, y, f(z, x)) + 2f(w, z, f(w, x)) - (6f(w, x, f(z, x))) - (2f(w, y, f(w, x))) - (6f(w, z, f(z, x)))
I’m only beginning to use SymbolicUtils.jl, so I don’t know if I’m using the right combination of rewriters.
This version is about 10 times slower than some Maple code I have (which is similar in spirit to expand_multilinear from @greatpet’s example above). Maybe something like expand_multilinear is the better approach, for example because the Fixpoint rewriter makes the Postwalk rewriter unnecessarily traverse the whole expression tree over and over again, instead of moving its way up from the bottom level only once. It would nevertheless be interesting to benchmark this Julia code also against @greatpet’s Mathematica example.
Regarding @acrule: From reading the documentation, I don’t fully understand what it does. It is not clear to me which functions are affected if I use @acrule for a pattern involving more than one function.
Regarding performance, I timed your code with a large-ish input,
r(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))
It took about 11 seconds on my laptop after a JIT warmup run. An equivalent Mathematica version runs more than 50 times faster, but interestingly the performance is similar to your Julia code if I run it with an open-source clone of Mathematica, expreduce (the software is very incomplete, but still good enough for this situation). I suppose the conclusion is just that Mathematica’s rewriting engine has been heavily optimized.
My experience with Maple is that due to a lack of built-in pattern matching capabilities, I have resort to custom code a lot, but unfortunately for loops are slow in Maple, so Julia is much nicer when you have to code up something from scratch.
OK, so this Julia code is too slow. I need to try something else. For the record, the following Maple code is at least 100 times faster than the Julia code on your example.
f := proc()
local a, i, x, y;
a := [args];
for i to nargs do
x := a[i];
if type(x, `+`)
then return `+`(seq(procname(op(1..i-1, a), y, op(i+1..-1, a)), y = x))
elif type(x, `*`)
then return `*`(op(1..-2, x), procname(op(1..i-1, a), op(-1, x), op(i+1..-1, a)))
fi
od;
'procname'(args)
end;
It would be great to have a Julia solution with a performance comparable to Maple and Mathematica.
Here’s a hand-coded Julia version which is 20 times faster than the version using rewrite rules, though it’ll need further optimization to be competitive. It’s a little long, but I’ll just dump the 60 lines of code here.
using SymbolicUtils
import SymbolicUtils.Add, SymbolicUtils.Mul, SymbolicUtils.Term
@syms u v w x y z argument_position(a,b)
f = SymbolicUtils.Sym{(SymbolicUtils.FnType){Tuple, Number}}(:f)
function expand1(term)
if !isa(term, Term) || !(term.f === f)
return term
end
args = deepcopy(term.arguments)
for i in 1:length(args)
if isa(args[i], Add)
added = summands(args[i])
args[i] = sum([argument_position(i, added[j]) for j in 1:length(added)])
else
args[i] = argument_position(i, args[i])
end
end
# For example, for term=f(x+2y, z), so far
# args == [argument_position(1,x)*argument_position(2,z) + argument_position(1,2y)*argument_position(2,z)
allterms = summands(expand(*(args...))) # summands() is defined below
sum(map(expand_multilinear, map(to_function, allterms))) # to_function() and expand_multilinear() are defined below
end
# For example, summands(x+2y) == [x, 2y]
function summands(x::Add)
d = x.dict
[term[1]*term[2] for term in d]
end
# Any other type
summands(x) = [x]
# For example, to_function(argument_position(1, x) * argument_position(2, y)) == f(x,y)
function to_function(x::Mul)
d = x.dict
args = Vector{Any}(undef, length(d))
for term in d
@assert term[2] == 1
@assert term[1].f == argument_position
args[term[1].arguments[1]] = term[1].arguments[2]
end
f(args...)
end
# For example, expand_multilinear(f(x,2y)) == 2f(x,y)
function expand_multilinear(expr)
function_arguments = Vector{Any}(deepcopy(expr.arguments))
overall_coeff = 1
for i in 1:length(function_arguments)
# check if this function argument is a product of things
if isa(function_arguments[i], Mul)
current_coeff = function_arguments[i].coeff
# absorb nontrivial coefficient into an overall factor
if current_coeff != 1
overall_coeff *= current_coeff
function_arguments[i] = function_arguments[i] / current_coeff
end
end
end
overall_coeff * f(function_arguments...)
end
@time expand1(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));
# It takes about 0.64 seconds on my laptop
Here is a version that only needs 115ms for your example, just 20% slower than Maple. (My laptop seems to be about as fast as yours.) So one can get a performance similar to Maple and Mathematica, and maybe even better with additional tweaking.
So far, the multilinear function f is hard-coded, and there may be other things one could improve. Note that all occurrences of f in an expression are expanded:
julia> expand_multilinear( f(f(x+2*y), f(3*z-w)) )
3f(f(x), f(z)) + 6f(f(y), f(z)) - f(f(x), f(w)) - (2f(f(y), f(w)))
Code (updated)
using SymbolicUtils
using SymbolicUtils: Add, Mul, Term, Symbolic
const f = SymbolicUtils.Sym{(SymbolicUtils.FnType){Tuple, Number}}(:f)
function expand_multilinear_arg(t::Add, i::Int)
sum(c * expand_multilinear_arg(x, i) for (x, c) in t.dict)
end
function expand_multilinear_arg(t::Mul, i::Int)
a, _ = first(t.dict)
# by construction, t has only one factor, and that has exponent 1
t.coeff * expand_multilinear_arg(a, i)
end
function expand_multilinear_arg(t::Term{T}, i::Int) where T
f = t.f
b = convert(Vector{Any}, t.arguments)
# this is needed for the assignment of one(T) to b[i] later one
# note that the type of FnType.arguments is Any anyway
x = b[i]
if x isa Add
d = Dict((b[i] = y; f(b...)) => c for (y, c) in x.dict)
c = x.coeff
if !iszero(c)
b[i] = one(T)
d[f(b...)] = c
end
Add(T, zero(T), d)
elseif x isa Mul
p = x.dict
b[i] = length(p) == 1 ? first(p).first : Mul(T, one(T), p)
Mul(T, x.coeff, Dict(f(b...) => one(T)))
else
t
end
end
function expand_multilinear_allargs(t)
# the function f is hard-coded here
if t isa Term && t.f == f
foldl(expand_multilinear_arg, 1:length(t.arguments); init = t)
else
t
end
end
expand_multilinear = Rewriters.Postwalk(expand_multilinear_allargs)
3 Likes
A minor bug fix for your code: the line
y = length(p) == 1 ? first(p) : Mul(Number, 1, p)
should probably be
y = length(p) == 1 ? first(p)[1] : Mul(Number, 1, p)
That’s needed for getting correct results for expand_multilinear(f(2x, y)), at least on my machine. Anyway, I’m really excited about doing symbolic calculations in Julia. No Python package would offer competitive performance for this kind of user-programmed symbolic manipulations, unless someone implements an underlying pattern matching / term rewriting engine in C++, or Julia 
2 Likes
Thanks! I’ve updated the code. I’ve thought a bit about making it even faster, but I didn’t see anything.
EDIT: On larger examples the Julia/SymbolicUtils code above is again much slower than Maple. So there is still a long way to go.
Yeah, we’ll follow up in the performance in an issue. Can you share the benchmarks to Maple and Mathematica in the issue as well?
1 Like
I’ve added Julia/Maple benchmarks to the GitHub issue.
@greatpet: Could you add Mathematica timings?
1 Like
I just added the Mathematica timings to your Github issue (my Github name is zengmao).
2 Likes
That’s all super helpful, thanks!