Struggling with defining a good factorization rule in SymbolicUtils.jl

General Usage
1

I started playing with Symbolics.jl but started needing more rewriting rules than what’s currently proposed. Taking a quick look at the rewrite docs, I don’t understand why some of the examples below don’t work.

using SymbolicUtils

@syms a::Real b::Real c::Real
factor = @acrule (~x*~y + ~x*~z) => ~x * (~y + ~z)

factor(a*b + a*c)           # a*(b+c)
factor(sin(a)*b + sin(a)*c) # nothing
factor(2b+2c)               # 2(b+c)
factor(2a*b + 2a*c)         # nothing
factor(a^2 + a*c)           # nothing

Any clue as to what’s happening here and how to fix it? (I am sure there are much better ways to write a factorization rule, including an arbitrary number of sum terms which I couldn’t figure out with segment variables).

2

Try defining two rules.

factor1 = @rule (~x * ~y + ~x * ~z) => ~x * (~y + ~z)
factor2 = @rule (~x * ~z + ~y * ~z) => (~x + ~y) * ~z
factor1(a*b + a*c) # a*(b + c)
factor2(a*b + a*c) # nothing
factor1(sin(a)*b + sin(a)*c) # nothing
factor2(sin(a)*b + sin(a)*c) # (b + c)*sin(a)

I think the root problem is that SymbolicsUtils.jl is not aware of the inherent commutativity of the multiplication operator *. The @acrule macro only imposes associativity-commutativity for the top-level operator, which is the + operator for the rule with LHS (~x * ~y + ~x * ~z). This seems to be a serious limitation though not strictly speaking a bug. It’ll be great if Symbolics devs can offer some insights here. Should I file a GitHub issue?

The example factor(2a*b + 2a*c) will not work because 2*a*b cannot be matched by ~x*~y, a well documented behavior (or perhaps limitation) of SymbolicsUtils.

1 Like
3

Thank you for your answer! My guess woud be then that SymbolicUtils rearranges sin(a)*b into b*sin(a), thus why only right factorization works.

It would probably be a good idea to file an issue, I can do it if you lack the time.

Speaking of this, the rules seem to work on the top-level expression but never on sub-expressions: for instance (taking your factor1 and factor2),

factor = Rewriters.RestartedChain([factor1, factor2])
factor(2*(a*b + a*c)) # 2*(a*b + a*c)

Is this also part of SymbolicUtils’ limitations? Is there a way to circumvent this?

4

To match sub-expressions, use Prewalk or Postwalk:

using SymbolicUtils
factor1 = @rule (~x * ~y + ~x * ~z) => ~x * (~y + ~z)
factor2 = @rule (~x * ~z + ~y * ~z) => (~x + ~y) * ~z
factor = Rewriters.Prewalk(Rewriters.Chain([factor1, factor2]))
factor(2*(a*b + a*c)) # 2a*(b + c)

Indeed I’m a bit too busy to file an issue, though I wouldn’t expect a quick resolution anyway, since implementing the requested feature seems to be a difficult job.

1 Like
5

Thank you for the answer! Prewalk is indeed what I was looking for.

I created the issue for the AC problem – if you see ways to improve it don’t hesitate, I am not used to writing issues yet!

1 Like
6

Your description of the issue looks clear enough!

1 Like
7

I have a pending pull request which would solve this problem. With the new code, you can pass an extra fullac option to the @rule macro to enable fully nested associative-commutative pattern matching.

# warning: example for in-development code
julia> using SymbolicUtils

julia> @syms a b c

julia> factor = @rule (~x*~y + ~x*~z) => ~x * (~y + ~z) fullac;

julia> factor(sin(a)*b + sin(a)*c)
(b + c)*sin(a)

julia> factor(a*b + a*c)
a*(b + c)

To factor 2a*b + 2a*c, you’ll need to match a factor which consists of more than one term, using segment patterns (in combination with my new code):

julia> factor2 = @rule (~~x*~y + ~~x*~z) => *(~~x..., ~y+~z) fullac;

julia> factor2(2a*b + 2a*c)
2a*(b + c)

Here ~~x has matched with the arguments [2, a] of the * operator, and the list is splat on the RHS to recover the product 2a.

8

For Symbolics.jl, it is possible to use the new factor_use_nemo (uses Nemo.jl):

julia> using Symbolics, Nemo

julia> @variables a b c;

julia> f = expand((a - 1)*(a^2 + b^2)*(a - c));

julia> Symbolics.factor_use_nemo(f)[2]
3-element Vector{Num}:
 -(1//1) + a
   a^2 + b^2
       a - c

We might want to export a version of it.

2 Likes
9

That’s very nice. I’d love to see an additional feature: it should be possible to use the Nemo backend to factorize more general expressions, e.g. log(a)^2 - log(b)^2 == (log(a) + log(b)) * (log(a) - log(b)), if you map the two logarithms (or any other functions) into polynomial variables.

I was just using the factorization example to demonstrate associative-commutative pattern matching. Of course a serious library like Nemo is the way to go to do factorizations properly.

1 Like