MatrixSymbol, exp, diff, and simplification

[saupin guillaume]

Hello,

I'm trying to use MatrixSymbol to do some computation.

More precisely, I have a function f(q1, q2) = e(A*q1) * e(B*q2), where

* q1, q2 are real

* A, B are 3x3 matrices

* e(M) is the matrix exponential

I want to compute  F = df(q1, q2) / dq2 * f(q1, q2)**-1 

F = e(A*q1) * B*e(B*q2) * (e(A*q1) * e(B*q2))**-1

F = e(A*q1) * B*e(B*q2) * (e(-B*q2) * e(-A*q1)

that can be simplified to 

F = e(A*q1) * B * e(-A*q1)

I tried to make these simple computations using sympy and the MatrixSymbol. I had to make some modification inside core/mul.py so that the diff() method works on MatrixSymbol.

With the followind code I almost manage to get a satisfying resuls, except for the final simplifiaction that does not work :(

from sympy import *

from sympy.matrices import *

q1 = Symbol('q1')

q2 = Symbol('q2')

q3 = Symbol('q3')

A = MatrixSymbol('A', 3, 3)

B = MatrixSymbol('B', 3, 3)

I = Identity(3)

print (A*B)**-1

print (exp(q1*A)*exp(q2*B))**-1 # the inversion is not computed

print (exp(q1*A)*exp(q2*B) * I)**-1 #inversion is computed if we add I

print (exp(q1*A)*exp(q2*B) * I).diff(q1) * (exp(q1*A)*exp(q2*B) * I)**-1 # The expression is not simplified-1

Any ideas ?

--
guillaume saupin

--
guillaume saupin

[Matthew Rocklin]

So, the simple answer is that SymPy doesn't currently support matrix exponentials. I.e. we don't support expressions of the type a**B where B is a matrix. 

When designing matrix expressions a choice was made to subclass from the normal expression class. This means that a lot of things from normal sympy work perfectly on matrix expressions without any additional code. Unfortunately some things that should raise errors go through, treating the matrix as a normal scalar symbol. Your case with exp(B) is such an example. 

I've created an issue for this

http://code.google.com/p/sympy/issues/detail?id=3119

For your immediate problem it's possible you could find a solution. You've gotten surprisingly far already. Whenever you multiply by I you're forcing sympy to turn the expression back into a matrix. Are you working off of a git branch somewhere? 

-Matt

--
guillaume saupin

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sympy?hl=en.

[kayhman]

On 29 fév, 19:53, Matthew Rocklin <mrock...@gmail.com> wrote:
> So, the simple answer is that SymPy doesn't currently support matrix
> exponentials. I.e. we don't support expressions of the type a**B where B is
> a matrix.

I've guessed that reading the code. Would it be hard for me to
implement the missing MatPow function ?

> <http://code.google.com/p/sympy/issues/detail?id=3119>For your immediate
> problem it's possible you could find a solution. You've gotten surprisingly
> far already. Whenever you multiply by I you're forcing sympy to turn the
> expression back into a matrix.

Would it be possible to create a specific exp function, that inherits
from ExpBase, or Exp, and that would return expr that have
MatrixSymbol type ?

It may avoid to add the "I" matrix to coerce expression into the
matrixSymbol type, no ? And hence the simplification should work ?

> Are you working off of a git branch
> somewhere?
>

Following your coding guidelines, I've created a github repository :
g...@github.com:kayhman/sympy.git.

[Matthew Rocklin]

I've guessed that reading the code. Would it be hard for me to
implement the missing MatPow function ?

 

It should be easy to implement such a class. I would suggest the name MatExponential. We already have a MatPow for the case A**b (A a matrix, b a scalar). I would suggest looking at MatPow and MatMul as examples. You'll also have to change the __rpow__ operator in MatExpr (look at the surrounding functions to see what to do). This should allow for the syntax e**A but not exp(A).

 

> <http://code.google.com/p/sympy/issues/detail?id=3119>For your immediate
> problem it's possible you could find a solution. You've gotten surprisingly
> far already. Whenever you multiply by I you're forcing sympy to turn the
> expression back into a matrix.

Would it be possible to create a specific exp function, that inherits
from ExpBase, or Exp, and that would return expr that have
MatrixSymbol type ?

 

You could make some other exp function like "mexp" or you could add a couple of lines in the exp function to check for matrix inputs. In general it is best to avoid changing core to help matrix expressions. I personally prefer the e**A method.

Following your coding guidelines, I've created a github repository :
g...@github.com:kayhman/sympy.git.

Excellent. The SymPy community uses git almost exclusively. This makes it much easier for us to see, experiment with, and comment on your code. 

[Ronan Lamy]

Le jeudi 01 mars 2012 à 08:20 -0600, Matthew Rocklin a écrit :
> I've guessed that reading the code. Would it be hard for me to
> implement the missing MatPow function ?
>
> It should be easy to implement such a class. I would suggest the name
> MatExponential. We already have a MatPow for the case A**b (A a
> matrix, b a scalar). I would suggest looking at MatPow and MatMul as
> examples. You'll also have to change the __rpow__ operator in MatExpr
> (look at the surrounding functions to see what to do). This should
> allow for the syntax e**A but not exp(A).

x**A is a confusing notation and is very rarely used with x != E. We
should only have exp(A).

[kayhman]

> x**A is a confusing notation and is very rarely used with x != E. We
> should only have exp(A).

I agree with that, but the solution proposed by Matthew seems more
simple to implement, at least for me.

Maybe when I get more accustomed to sympy I'll try to implement this
solution.

[kayhman]

> It should be easy to implement such a class. I would suggest the name
> MatExponential. We already have a MatPow for the case A**b (A a matrix, b a
> scalar). I would suggest looking at MatPow and MatMul as examples. You'll
> also have to change the __rpow__ operator in MatExpr (look at the
> surrounding functions to see what to do). This should allow for the syntax
> e**A but not exp(A).

I've implemented such a class, and it seems to works. I can now manage
to compute (e**(q*A)).diff(q).

However, I still have a major problem : sympy does not manage to
simplify e**A * (e**A)**-1, though it can simplify e**q * (e**q)**-1
where A is a SymbolMatrix, and q a symbol.

See the code below. Why is that ?

from sympy import *
from sympy.matrices import *

e = Symbol('e')

q1 = Symbol('q1')
q2 = Symbol('q2')
q3 = Symbol('q3')


A = MatrixSymbol('A', 3, 3)
B = MatrixSymbol('B', 3, 3)

print (e**(q2*B)) * (e**(q2*B))**-1 # not simplified :(
print (e**(q2*q1)) * (e**(q2*q1))**-1 # simplified to 1

[Aaron Meurer]

I don't see any reason to not support it. Maybe we could
automatically convert it to exp(log(x)*A) in this case, to make it
clear, though personally I don't even see why we would have to do
that. If no one uses x != E, then it won't ever be used.

Aaron Meurer

[Aaron Meurer]

Out of curiosity, why do you use Symbol('e') instead of E (i.e.,
exp(1))? They both behave the same here, but I think the latter is
what you really want.

I think it might have something to do with how MatrixExprs are not
automatically simplified in the core beyond those things the core
already knows about. Matthew will have to comment more. The solution
might be to make exp(A)**-1 automatically convert into exp(-A) like it
does in the normal core. I tried that and it seems there's a bug in
your branch:

In [22]: (e**(q2*B)) * (e**(-q2*B))
Out[22]: e

I also noticed that there's some strange pretty printing with this:

In [11]: (e**(q2*B)) * (e**(q2*B))**-1
Out[11]:
q₂⋅B ⎛ q₂⋅B⎞
e ⋅⎝e ⎠^-1

And by the way, to make it so you can use exp() instead of E**, I
think the best way would be to do what you're doing and fix
http://code.google.com/p/sympy/issues/detail?id=1799.

Aaron Meurer