Manual control of simplification and evaluation

[jtrakk]

Currently SymbolicUtils uses constructor-level simplification and automatic evaluation of expressions. This may be convenient for numerical work, but makes SymbolicUtils inappropriate for purely symbolic manipulation. It also eliminates algorithmic speedup opportunities.

There is an asymmetry: a raw symbolic expression can always be simplified for numerical work, but once simplified an expression cannot be unsimplified for symbolic work. To support both kinds of users and ensure high performance, I would like to propose that SymbolicUtils adopt a posture of not simplifying or evaluating without explicit user request.

For example, with @syms a b:

• a + 1 would stay a + 1, not changing to 1 + a as it currently does.
• substitute(a+b, Dict(a=>1, b=>2)) would go to Term(+, [1,2]), not 3 as it currently does.

The current SymPy maintainer has warned of this issue's critical importance. He says [emphasis added]:

Personally I think that automatic evaluation is the major mistake made in sympy (and many other CAS). That is what currently slows down ordinary operations in sympy. Automatic evaluation is extremely hard to change when backwards compatibility comes into play because any change is not compatible.

I hope we can get these decisions built into the early foundation of Symbolics.jl so we don't follow the wrong path.

[ChrisRackauckas]

a + 1 would stay a + 1, not changing to 1 + a as it currently does.
substitute(a+b, Dict(a=>1, b=>2)) would go to Term(+, [1,2]), not 3 as it currently does.

No, the data structures are different in SymbolicUtils.jl. The no-evaluation form is the simplified form in this case, by design, and that's why it's fast. There is no call to simplify required to make that stuff work, and it would be slower to hold it as a tree instead of 3. So this seems to be just a major misunderstanding of the design. This is a similar design to symengine.

[jtrakk]

That automatic simplification is exactly the design that I'm concerned about. Quoting another CAS developer in that thread

Simple, unevaluated expressions by default are a win-win: easier to implement, easier and more versatile to use (the user can decide precisely which transformations to apply), and more efficient (symbolic expressions are a terrible data structure for algebraic computation). Instead of performing clever and invariably buggy automatic “canonicalization” on expressions, a far better solution is to convert to specialized data structures internally for computations (e.g. polynomial types for expanding or factoring).

Instead of automatically canonicalizing a + 1 to 1 + a, it could just stay how it is until the user requests a transformation.

[ChrisRackauckas]

It's not canonicalized to either.

[jtrakk]

julia> using SymbolicUtils # v0.13.3

julia> @syms a
(a,)

julia> a + 1
1 + a

Currently a + 1 swaps its operands to produce 1 + a (unless this is just a quirk of the repr?). My proposal is that it not do that and to retain the input a + 1.

[ChrisRackauckas]

No it does not. That's just the printing. Look at the data structure. It is neither direction because it's not ordered.

[jtrakk]

julia> a + 1
1 + a

help?> SymbolicUtils.Add
  Add(T, coeff, dict::Dict)


  Represents coeff + (key1 * val1) + (key2 * val2) + ...

  where keys and values come from the dictionary (dict). where coeff and the vals are <:Number and keys are symbolic.

    •  operation(::Add) – returns +.

    •  symtype(::Add) – returns T.

    •  arguments(::Add) – returns a totally ordered vector of arguments. i.e. [coeff, keyM*valM, keyN*valN...]

The canonicalization of a + 1 in this case eliminates the order of the arguments. This transformation may be helpful in many cases, but I'd like to decide when to apply it. Likewise, I would like a purely symbolic substitution: substitute(a+b, Dict(a=>1, b=>2)) to 1+2, not to 3.

The symbolic semantics I'm describing are useful for

• manipulating symbols without making assumptions about the meaning of those symbols
• performance in symbolic computation

CCing some other people who may have thoughts @0x0f0f0f @fredrik-johansson @JacquesCarette @oscarbenjamin.

It might be helpful to have examples of the performance benefit if somebody knows a good one off hand.

[shashi]

There is an obvious tension between taking expressions as-is and having efficient datastructures which can make some assumptions.

For variables of numeric type (by default in @syms a -- a is of the symbolic type Number), we use the datastructure which makes it faster, smaller by utilizing the commutativity and associativity properties. Arguably it would not work for * and Quaternions but in that case, you can special case it.

The benefit of the multiple orders of magnitude of speed gained by this cannot be understated -- it makes new things possible.

It seems like there is room for a mode where you can treat expressions as is. We easily add a macro for this.

@macroexpand @symquote  a + 1

:(term(+, a, 1))

[shashi]

substitute(a+b, Dict(a=>1, b=>2)) to 1+2, not to 3.

This is fixed now (will be releasing), if you pass in fold=false.

[ChrisRackauckas]

It might be helpful to have examples of the performance benefit if somebody knows a good one off hand.

You can check the old PRs on this which documented the performance. It was a good 1000x speedup or so on the examples we had on-hand, and similar speedups were seen on biopharma real-world examples. At this point it's pretty much the other way around: there is so much evidence about how much of an acceleration that it gives in real examples that we would need to see some very compelling evidence to change to another form. I was actually pretty reluctant about it at first, but the results don't lie. So try something on say the BCR reaction from ReactionNetworkImporters, or on the robot mass matrix example, and see if you can find another form.

[YingboMa]

Do you have any concrete examples of when canonicalization is bad? It's not like change a + 1 to 1 + a does any expensive computation nor does it hinder any transformation of symbolic expressions. On the contrary, canonicalization massive speed up many passes.

[shashi]

You can also try @syms x::Any this will give you a blank-slate symbolic variable, and of course + etc are not defined on it.

[jtrakk]

@shashi Thanks for recognizing the representation tension.

substitute(a+b, Dict(a=>1, b=>2)) to 1+2, not to 3.
This is fixed now (will be releasing), if you pass in fold=false.

On master (efaa595) right now it returns 3:

julia> substitute(a+b, Dict(a=>1,b=>2), fold=false)
3

This function is more like an evaluate than a CAS subsitute.

You can also try @syms x::Any this will give you a blank-slate symbolic variable, and of course + etc are not defined on it.

julia> @syms a::Any b::Any
(a, b)

julia> a + b
ERROR: MethodError: no method matching +(::SymbolicUtils.Sym{Any, Nothing}, ::SymbolicUtils.Sym{Any, Nothing})

That's not what I have in mind: the result I want is just a + b.

The @symquote behavior you mention would be cool, but it would be useful as the default behavior, for functions I don't control.

@YingboMa

Do you have any concrete examples of when canonicalization is bad? It's not like change a + 1 to 1 + a does any expensive computation nor does it hinder any transformation of symbolic expressions. On the contrary, canonicalization massive speed up many passes.

I'll leave that to the CAS experts I mentioned.

[oscarbenjamin]

Do you have any concrete examples of when canonicalization is bad?

For some applications you really want to have control over how an expression is manipulated. Symbolic manipulation has applications that are not just about computing things and in which fine control is needed. For example you might use a CAS to run through the steps in a complicated derivation that is used in a mathematical text. Then if your expressions can be converted to latex/mathml/etc you can have the equations of a document generated semi-automatically. In this context you really don't want canonicalisation to mess up your equations though. If canonicalisation is unavoidably built in to core operations like substitution then the CAS is unusable for this.

Canonicalisation is also very much a slippery slope. If implicit/automatic canonicalisation is expected by users then there will be a long tail of different kinds of canonicalisation that some of them will want e.g. 1/sqrt(2) -> sqrt(2)/2 and so on. For every canonicalisation that one person wants there will be others who really don't want it because it then becomes impossible to represent expressions in different forms without auto-canonicalising.

It's not like change a + 1 to 1 + a does any expensive computation nor does it hinder any transformation of symbolic expressions.

If you want to think about expensive computation then you need to think about large expressions and/or operations that can happen many times. If your sum has millions of terms do they need to be sorted every time any substitution is performed? What if the terms themselves are complicated expressions and comparing two terms under the sorting order is nontrivial?

On the contrary, canonicalization massive speed up many passes.

And it can also waste time for no reason in situations that don't need it. What if you have an algorithm that works by repeatedly making many different substitutions? Perhaps a single canonicalisation at the end would have been fine but instead the sort has to be recomputed at every step because canonicalisation has been made automatic and unavoidable.

Ultimately high-level routines built on top of primitive routines like substitute can be made more efficient if they can control how operations at the lower level are handled. Every time something like canonicalisation is made automatic the higher level routines lose that control.

I would advocate for isolating primitives as much as possible so e.g. substitute should just perform the substitution and nothing more. If the output should be canonicalised then a separate canonicalise function should be used. Of course there can be convenience functions that combine them but there needs to be clear primitives that don't mix operations together. Having clear primitives gives control to users and makes a good basis for efficient higher-level operations.

[ChrisRackauckas]

This discussion does not seem to be relevant to the actual data structures of the repository. Could you give a real-world example, say an SBML file, which is handled better with a different data structure?

[JacquesCarette]

This is a very complex problem. It's not that dis-similar to choosing between call-by-value or call-by-name (or call-by-need); there are examples of each case that are clear win/lose (depending on what you're trying to prove). [And @oscarbenjamin 's answer came in while I was writing... I agree with all he says. So I'll assume that, and comment further.]

There is definite tension between naive usability and large-scale computation. Having + take care of both is where things kind of fall down. You want to be able to represent addition in a way where there are 'no surprises' (i.e. full control), and another where some amount of evaluation can happen. My current take is to provide 3 operations:

1. a fully inert + that does nothing at all, except be a representation of a sum
2. a + that does a variety of simplifications (which ones should be controllable with options)
3. a user-level + that's built on top of the first two, carefully documented on how it is built.

The default would be to use the user-level + until more control is definitely needed. The big mistake is to conflate all 3 into a single command.

[oscarbenjamin]

Could you give a real-world example, say an SBML file, which is handled better with a different data structure?

I'm not personally experienced enough with Julia to do this (I don't even know what a SBML file is).

I've been trying Symbolics/SymbolicUtils out but so far I'm hitting against quite basic things like JuliaSymbolics/Symbolics.jl#328 (comment)

[shashi]

On master (efaa595) right now it returns 3:

julia> substitute(a+b, Dict(a=>1,b=>2), fold=false)
3

Oh right... This is because of how similarterm works on the closure of types that we use to represent expressions of type <:Number by default. JuliaSymbolics/Symbolics.jl#328 (comment) should work. But the inconsistency is not lost on me. My proposed "proper" solution would be as below:

In Julia there are not multiple + generic functions. The right way to implement an as-is symbolics would be to define a different closure of symbolic types. This could be a module,

AsIsModule.@syms a b c

could give you a different type of symbol. Then the + generic function can be defined on that to return a term type that does not canonicalize. And then you'd need to implement similarterm on these two types to keep the answer within these two types.
This should be easy to implement (I would say 200-300 lines of code including the macro).
Then substitute should actually work as you describe without any changes to substitute's code. This is a larger goal of SymbolicUtils: one should be able to decouple symbolic manipulation code and the representation and how the default operators work on those representations.

If you go through the history of this package you will notice that we did start off with this. But we found that it is really awful for systems with 30000 equations in 10000 variables, and generating code for that. It was unreasonable to wait for them to simplify. Switching to the canonical form lead to saving hours.

I have some thoughts about the social aspects and organizing symbolic programming to support all use cases. I may write a post about this. But I would like to quickly note here: In other ecosystems, we are conditioned to think "one package -- one representation" and any mis-step with it is decades of lock-in. But reality is much more forgiving when you are in a language with dynamic multiple dispatch.

[shashi]

Closed by #429