simplify the expression after multiplying matrix of complex numbers

[Ilshat Garipov]

I recently faced with the need to do small calculations with matrices opf complex number in any symbolic math environment, so I decided to try with simpy. And here is an example of my code:

import sympy as sym
from sympy import *
from sympy.physics.quantum.dagger import Dagger

a1, a2, a3, a4 = sym.symbols('a1, a2, a3, a4', real=True)
b1, b2, b3, b4 = sym.symbols('b1, b2, b3, b4', real=True)
g1, g2, g3 = sym.symbols('g1, g2, g3', real=True)
h1, h2, h3 = sym.symbols('h1, h2, h3', real=True)

U = sym.Matrix([[a1 + b1*1j, a2 + b2*1j], [a3 + b3*1j, a4 + b4*1j]])
G = sym.Matrix([[g1, g2 + g3*1j], [g2 - g3*1j, -g1]])
H = sym.Matrix([[h1, h2 + h3*1j], [h2 - h3*1j, -h1]])
U1 = Dagger(U)

Rg= U*G*U1
Rh = U*H*U1
tr = (Rg*Rh).trace()

The problem is that result looks awesomely bulky:

((a1 - 1.0*I*b1)*(g1*(a1 + 1.0*I*b1) + (a2 + 1.0*I*b2)*(g2 - 1.0*I*g3)) + (a2 - 1.0*I*b2)*(-g1*(a2 + 1.0*I*b2) + (a1 + 1.0*I*b1)*(g2 + 1.0*I*g3)))*((a1 - 1.0*I*b1)*(h1*(a1 + 1.0*I*b1) + (a2 + 1.0*I*b2)*(h2 - 1.0*I*h3)) + (a2 - 1.0*I*b2)*(-h1*(a2 + 1.0*I*b2) + (a1 + 1.0*I*b1)*(h2 + 1.0*I*h3))) + ((a1 - 1.0*I*b1)*(g1*(a3 + 1.0*I*b3) + (a4 + 1.0*I*b4)*(g2 - 1.0*I*g3)) + (a2 - 1.0*I*b2)*(-g1*(a4 + 1.0*I*b4) + (a3 + 1.0*I*b3)*(g2 + 1.0*I*g3)))*((a3 - 1.0*I*b3)*(h1*(a1 + 1.0*I*b1) + (a2 + 1.0*I*b2)*(h2 - 1.0*I*h3)) + (a4 - 1.0*I*b4)*(-h1*(a2 + 1.0*I*b2) + (a1 + 1.0*I*b1)*(h2 + 1.0*I*h3))) + ((a1 - 1.0*I*b1)*(h1*(a3 + 1.0*I*b3) + (a4 + 1.0*I*b4)*(h2 - 1.0*I*h3)) + (a2 - 1.0*I*b2)*(-h1*(a4 + 1.0*I*b4) + (a3 + 1.0*I*b3)*(h2 + 1.0*I*h3)))*((a3 - 1.0*I*b3)*(g1*(a1 + 1.0*I*b1) + (a2 + 1.0*I*b2)*(g2 - 1.0*I*g3)) + (a4 - 1.0*I*b4)*(-g1*(a2 + 1.0*I*b2) + (a1 + 1.0*I*b1)*(g2 + 1.0*I*g3))) + ((a3 - 1.0*I*b3)*(g1*(a3 + 1.0*I*b3) + (a4 + 1.0*I*b4)*(g2 - 1.0*I*g3)) + (a4 - 1.0*I*b4)*(-g1*(a4 + 1.0*I*b4) + (a3 + 1.0*I*b3)*(g2 + 1.0*I*g3)))*((a3 - 1.0*I*b3)*(h1*(a3 + 1.0*I*b3) + (a4 + 1.0*I*b4)*(h2 - 1.0*I*h3)) + (a4 - 1.0*I*b4)*(-h1*(a4 + 1.0*I*b4) + (a3 + 1.0*I*b3)*(h2 + 1.0*I*h3)))

And i am wondering if there are any built-in methods to shorten this expression? Am I doing any thing wrong in my code? I am new in python so there might be things I don't yet know.. Thanks in advance for helping.

[Oscar Benjamin]

It's better to use sympy.I rather than 1j as it is exact and does not
use floating point.

You say that the result is bulky but do you have a reason to think
that it can be any simpler?

This kind of explosion is typical if you have matrices where every
element is an arbitrary symbol. That's why it's useful to write
equations where each symbol is a whole matrix.

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/b6b78b13-f11d-427e-9633-26381d0f34fa%40googlegroups.com.

[Aaron Meurer]

The expression is basically already as compact as it can be. You can
use expand() to expand it into a single polynomial. The expression
itself is about six times bigger if you do that, but it will be
flattened out. You can see if you call factor() on the expression that
it doesn't factor (except for a common constant term of 8).

Aaron Meurer

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxSWbku%3D3Yx5j1%3DxfYON6WTAFvWmByXJj%2BA9Y7eJ2Me0TA%40mail.gmail.com.

[Chris Smith]

You can get a compact representation since there are many repeated expressions (though not appearing in a simple factorable manner):

>>> tr.count_ops()

271

>>> c1=cse(tr)

>>> count_ops(c1)

74

>>> c1[1]  # the expression with subexpressions extracted

[(x11*x24 + x23*x8)*(x2*x25 + x26*x9) + (x11*x9 + x2*x8)*(x14*x2 + x16*x9) + (x14*x23 + x16*x24)*(x2*x21 + x22*x9) + (x21*x23 + x22*x24)*(x23*x25 + x24*x26)]

[Ilshat Garipov]

  Hello Chris, thank you much for the tip! And thanks to all who responded to my issue. Very much appreciate it.   

I started to look for more information about count_ops()   and found the function - separatevars() - on the same page that did the trick . 

What I am trying to check is that if we have two special type complex matrix X and Y then trace of X*Y defines the Euclidean metric. 

And if R is  a mapping that translates each such matrix into a matrix U*X*U1 - then R is an ortogonal transformation, i.e if does preserve the metric. 

I initially made a mistake - not all Hermitian matrices will satisfy this property. But only unitary with determinant equal to 1. They have a common look like in a following solution:

import sympy as sym

from sympy import *
from sympy.physics.quantum.dagger import Dagger

z1, z2, w1, w2 = sym.symbols('z1, z2, w1, w2', real=True)
z = z1 + z2 * I
w = w1 + w2 * I
U = sym.Matrix([[z, w], [-conjugate(w), conjugate(z)]])
# https://docs.sympy.org/latest/modules/physics/quantum/dagger.html
U1 = Dagger(U)

x1, x2, x3 = sym.symbols('x1, x2, x3', real=True)
y1, y2, y3 = sym.symbols('y1, y2, y3', real=True)
X = sym.Matrix([[x1, x2 + x3 * I], [x2 - x3 * I, -x1]])
Y = sym.Matrix([[y1, y2 + y3 * I], [y2 - y3 * I, -y1]])

Rx = U * X * U1
Ry = U * Y * U1
F = Rx * Ry
tr = (Rx * Ry).trace()

print(separatevars(tr))

It gave me the following output (!) wow :  2*(x1*y1 + x2*y2 + x3*y3)*(w1**2 + w2**2 + z1**2 + z2**2)**2   (1)

unitary U with determinant equal to 1 means that w1**2 + w2**2 + z1**2 + z2**2 = 1.  (2)

So this means that result is 

2*(x1*y1 + x2*y2 + x3*y3)  (3)

which coincides with the metric = (X*Y).trace()


By the way - may anyone help with this -  if there is a  way to simplify an expression with respect to provided equation? 

For example, As in above - simplify (1) with respect to (2) gives (3) ?

Thank you!

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/cb2e2f76-22e3-46a4-af8c-274bc598ac00%40googlegroups.com.

[Oscar Benjamin]

Now that your result has a simple form you can also use factor rather
than separatevars.

For simplifying with respect to an equation you can use subs:

In [85]: factor(tr).subs(w1**2 + w2**2 + z1**2 + z2**2, 1)
Out[85]: 2⋅x₁⋅y₁ + 2⋅x₂⋅y₂ + 2⋅x₃⋅y₃

It is generally more reliable to substitute for as atomic a term as
possible though:

In [86]: factor(tr).subs(w1**2, 1 - w2**2 - z1**2 - z2**2)
Out[86]: 2⋅x₁⋅y₁ + 2⋅x₂⋅y₂ + 2⋅x₃⋅y₃

In [87]: factor(tr).subs(w1, sqrt(1 - w2**2 - z1**2 - z2**2))
Out[87]: 2⋅x₁⋅y₁ + 2⋅x₂⋅y₂ + 2⋅x₃⋅y₃

For simplifying rational functions with respect to polynomials there
is ratsimpmodprime that can take the equation directly as an expr:

In [90]: ratsimpmodprime(tr, [w1**2 + w2**2 + z1**2 + z2**2 - 1])
Out[90]: 2⋅x₁⋅y₁ + 2⋅x₂⋅y₂ + 2⋅x₃⋅y₃

(Note that factor does not need to be called first)

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAPTMStTx63stRK9inm4z0SktYVgKjfW-1DqTc0spB4terkDaeA%40mail.gmail.com.