Canonicalization - sympy/sympy GitHub Wiki

Canonicalization

This document will serve as a place to discuss ideas on how to better do canonicalization. In other words, the things that are currently done by things like Mul.flatten and Add.flatten. Feel free to edit it with your own ideas.

See also issue 1941.

Objects requiring custom logic

I think it would be helpful to list all the kinds of objects we can think of that require custom logic. Since Mul seems to be the most common case, I will just use it exclusively here.

Objects already handled in Mul.flatten

  • Infinities: oo, -oo, zoo, and nan. These objects generally eat up other objects, where it is mathematically correct. For example, oo*positive = oo, zoo*non-zero = zoo, nan*anything=nan.

  • Imaginary unit: I should combine with I to make -1 automatically. Integer powers of I are automatically evaluated, and complex numbers with real and imaginary parts are not automatically evaluated in a product, so I think this is the only case.

  • Non-commutatives: This is a tricky one, and in some sense, completely different from the rest. The only issue here is that the order should be preserved. Also, any commutative expression should be pull out. The way this is done right now is to keep track separately of the "commutative part" and the "non-commutative part". The commutative part is treated almost as a separate Mul multiplying the non-commutative part.

  • Numbers: Specifically, Integers, Rationals, and Floats. These are always combined into a single number. In Mul.flatten, they are treated separately as the coefficient.

  • Rational powers of rational numbers: Right now, we apply a canonicalization algorithm on rational powers of rational numbers, which reduces powers of the factors of the numbers (we don't factor very large integers, but this is irrelevant here I think). For example, (63/16)**(5/6), 3**(5/3)*(7/16)**(5/6), and (63)**(5/6)/(8*2**(1/3)) are all changed into 3*6**(2/3)*7**(5/6)/16. Once again, all terms must be gathered for this to happen correctly. For example, sqrt(2) and sqrt(6) by themselves to not simplify, but sqrt(2)*sqrt(6) simplfies to 2*sqrt(3). See for example issue 415, and the fix that had to be applied there.

  • Order Terms: Order terms (like O(x**2)) are currently special cased in both Mul.flatten and Add.flatten so that things automatically combine with them. So, e.g., O(x**2)*x becomes O(x**3), O(x**2) + x**3 becomes just O(x**2).

Objects not handled in Mul.flatten, but that would like to be

  • Arbitrary constants: See issue 1336. It would be nice to have a special type of Symbol that automatically absorbs constants into it. This could be (optionally) returned from dsolve() and integrate(). In dsolve(), this would eliminate the need for constantsimp(). There are already many tests for this in sympy/solvers/tests/test_constantsimp.py. This is very similar to how O() works.

  • A ComplexNumber type, which would act just like Number + Number*I except ComplexNumber*ComplexNumber would automatically expand and give another ComplexNumber.

  • MatrixExpr: A MatrixExpr should be only able to multiply another MatrixExpr if their inner shapes match. They should only be able to add if their shapes match. Matrix multiplication is non-commutative. Any regular Expr (or at least a commutative one) should be able to multiply a MatrixExpr and be pulled out as the coefficient.

  • Block Matrices:

  • ImmutibleMatrix:

  • Quantum Stuff:

  • Units and quantities:

  • Differential Operators: Using multiplicative notation. For example, p = DerivativeOperator(x), p*sin(x) => cos(x).

Problems

  • Commutativity means that we have to consider how an object combines with everything in the Mul, not just what it is next to. So, for example, zoo*x does not reduce to zoo by default, because x could be 0 (in which case it would become nan). But zoo*x*3 should reduce to zoo*x, i.e., the 3 should absorb into the zoo. How can we use Python's double dispatch along with this?

  • Commutativity also means that the order of the final expression can (and should) be canonicalized, so that it's easy to assert that x*y == y*x.

  • Associativity holds for anything in a Mul, but if rules conflict with each other x*(y*z) might not give the same thing as (x*y)*z, because of the order that things are evaluated in by Python (or by whatever method we end up using). Is this a problem even with non-conflicting rules? How should we deal with conflicting rules?

Unevaluated expressions

  • What degree of non-evaluation should we allow? Should we allow 1*1*1*1 (currently impossible). Should we allow to create x*y and y*x as different objects? What about x*(y*z) and (x*y)*z? How are these expressions treated by other functions, which implicitly expect canonicalization to have occurred.

Resuming one possible solution taken from mail-list discussions

See: https://groups.google.com/d/topic/sympy/fCQEdSQybTM/discussion

See: https://groups.google.com/d/msg/sympy/rKnqkU_iK44/lzMBjSfN6EcJ

See: https://groups.google.com/d/topic/sympy/kdhLn36ZpRE/discussion

Double dispatch rules used only for operator precedence, not for simplification. That means that x+x+y*z is evaluated to Add(x,Add(x,Mul(y,z))). No associativity or commutativity or anything. After the AST is created a default canonicalizer is called (i.e. all the logic from .flatten()). This way one needs not subclass Expr (e.g. MatrixSymbol or the quantum module). Also instead of writing algorithms for each operation one can write rules executed by the same algorithm.

Graph Optimizers - Theano's solution

Theano is an open source project that serves as a bridge between symbolic and numeric computation. Like us, they build up an symbolic expression and simplify it. Unlike us they have clearly separated these two processes. They have a nice way of specifying individual simplification rules and applying them sequentially. Looking into their system may provide us with some ideas.

Some links: Theano main page, optimizers, graph optimization.