Associative noncommutative rules

[david-pl]

I'm working with q-commutative algebra, similar to #96, where some variables are noncommutative under multiplication but follow fixed commutation relations. In order to specify these, it would be practical to have associative noncommutative rules. These would be similar to ACRule but keep the order of the arguments.

For example, consider two noncommutative variables a::TypeA and b::TypeB, which have the commutation relation [a,b]=1. I would like to rewrite a*b => b*a +1. If I do

@rule *(~x::sym_isa(TypeA), ~y::sym_isa(TypeB)) => ~y*~x + one(~x)

then a*b is simplified correctly, but as soon as there are multiple arguments such as in 2*a*b it does no longer work. What I currently do is step through the arguments of a multiplication and look for consecutive occurrences of a and b.
This should be much simpler with something like

@ancrule *(~x::sym_isa(TypeA), ~y::sym_isa(TypeB)) => ~y*~x + one(~x)

I found that this can be almost done already with the following:

subsets(inds, n) = (@views(inds[i:i+n-1]) for i=1:length(inds)-n+1)
macro ancrule(expr)
    arity = length(expr.args[2].args[2:end])
    quote
        ACRule(subsets, $(esc(:(@rule($(expr))))), $arity)
    end
end

The only problem with this approach is that when simplifying, this line

SymbolicUtils.jl/src/rule.jl

Line 326 in 3fd27d8

return @timer "acrule" Term{T}(f, [result, (args[i] for i in eachindex(args) if i ∉ inds)...])

does not keep the order in the result. This leads to issues such as infinite recursion in expressions like b*a*b*a. Changing that to keep the order works for my case, but causes other issues (tests fail).

Alternatively, one can implement a custom ANCRule type, which is essentially a copy-paste of ACRule, but with a slightly modified application to expressions. This is what I found to work best for now, since it also doesn't require changes in SymbolicUtils directly.

Would it make sense to have this implemented directly in SymbolicUtils? Is this even the right approach? Any input is appreciated, thanks.

[shashi]

Nice issue! Thanks for digging for the problem!

Maybe we could store another callback in ACRule which is called with the result and match and returns the final term? This comprehension could be what the default callback does.

[shashi]

Do you think it might work better if you just used @rule?

@rule *(~~a, ~x::sym_isa(TypeA), ~y::sym_isa(TypeB), ~~b) => *((~~a)..., ~y*~x + one(~x), (~~b)...)

Yup that should do the right thing!

[david-pl]

@rule *(~~a, ~x::sym_isa(TypeA), ~y::sym_isa(TypeB), ~~b) => *((~~a)..., ~y*~x + one(~x), (~~b)...)

Thanks, this works perfectly! Still need some practice in finding what the "best" way to write things down is. Thanks a lot!

[shashi]

Haha the rule of thumb is @rule should be able to do anything that the other rules can do (may require multiple @rules, but will work).

I'd welcome any documentation on this! Maybe a separate markdown file in docs/ with some thoughts on how to write rules for Q-Algebraic operations like these. Let me know if you end up using this in a package. Would love to have some examples we can link to!

[david-pl]

I was able to do it with @rule before, but I stepped through all the arguments of an expression *(~~a) and looked for consecutive occurrences. So in a sense I wrote my own matcher, which is silly of course. It worked, but your way is much more elegant.

Well, I am actually writing a package using this. It's not released yet, but hopefully will be at some point. It's in the context of quantum optics. The basic idea is to derive equations of motion for Q-commutative variables. Then you can perform some approximation (cumulant expansion) to convert these equations to c-number equations. Finally, I generate a function from this (basically just like MTK), which can be used in OrdinaryDiffEq.
Here's the link: https://github.com/david-pl/Qumulants.jl

Most of the rules there are just copy-paste from SymbolicUtils, actually. You only need your own sorting (keeping non-commutativity in mind), and you need to expand + inside of * in order to match pairwise commutation relations. I'll see if I can write a nice MWE illustrating this for the documentation.

[ChrisRackauckas]

Out of curiosity, is there a reason Qumulants isn't using MTK? Is it because of MTK can't work with the alternative algebras? This looks like something where you'll want to have features like sparse multithreaded Jacobians, so I hope we don't end up with diverging implementations of all of that, so I'd be curious to figure out where the abstractions have gone wrong (and it might just be in the algebras portion)

[david-pl]

It's just because of the non-commutative algebra Qumulants needs. At some point I'll switch to MTK for generating ODE functions. Right now this part of Qumulants is in a very basic state - I just wanted a cheap implementation that works for now. I don't plan on re-implementing things that MTK does really well anyways. Taking the current state of Qumulants, I could already convert the c-number expressions to MTK and leave things such as code generation and computing Jacobians to it.

Anyway, down the line Qumulants will ideally just implement the symbolic derivation of Heisenberg equations (non-commutative differential equations) using a set of commutator rules, and functionality to convert these into commutative ones. Everything else would be SymbolicUtils and MTK.

[ChrisRackauckas]

Sounds like a great application. I'll follow what you do and hopefully we can get you to where you need to be. Feel free to liberally file issues in both libraries.

[shashi]

You can call ModelingToolkit.to_mtk(expr) to turn a SymbolicUtils expression into ModelingToolkit one. You may have to declare function-like symbols though.

@variables t a b(t)

will need to be:

@syms  t() a() b(t)
a = a()
b = b(t())