Matrix Simplification
Addison Cugini
I'm working on implementing a simulation of a quantum computer in Sympy. One thing I have to do is make a matrix representation of a series of quantum gates and am having some troubles with the simplifications that automatically occur. I basically have two things that result in hard to read matrix expressions:
1. Sympy automatically turns exp(I*pi/2) into i and exp(I*pi) into -1, which results in ugly expressions like I*exp(I*pi/4) rather than the cleaner expression exp(3*I*pi/4).
2. I can't find a way to factor out constants that are common to all elements in a matrix. (In fact I can't even create a Mul object which holds a matrix and a constant because a sympy-matrix cannot be simpified)
I've attached a screen-shot of one of the rather unpleasant looking matrices. Is there any way to fix either one of these problems?
Thanks,
Addison J Cugini
Aaron S. Meurer
I think Chris had a branch somewhere with some kind of extract_gcd() method for Matrix. Maybe it could help you.
By the way, Maple automatically reduces exp(3*I*pi/4) to -sqrt(2)/2 + I*sqrt(2)/2. So maybe SymPy's exp() needs to be smarter so it can completely reduce itself, instead of only partially. Ideally, I think, it would use roots() or some equivalent to find the relevant roots of unity and replace it with those if they are radical.
Aaron Meurer
Brian Granger
On Mon, Jul 12, 2010 at 12:26 PM, Aaron S. Meurer <asme...@gmail.com> wrote:
> Well, fixing sympify to work with matrices would be a good start. There would probably be more work needed than that to make Matrix play nice with Mul, though.
If you end up tackling this, make sure you do this work in a new
branch from a fully updated master branch.
Cheers,
Brian
> I think Chris had a branch somewhere with some kind of extract_gcd() method for Matrix. Maybe it could help you.
>
> By the way, Maple automatically reduces exp(3*I*pi/4) to -sqrt(2)/2 + I*sqrt(2)/2. So maybe SymPy's exp() needs to be smarter so it can completely reduce itself, instead of only partially. Ideally, I think, it would use roots() or some equivalent to find the relevant roots of unity and replace it with those if they are radical.
>
> Aaron Meurer
> On Jul 12, 2010, at 12:41 PM, Addison Cugini wrote:
>
>> Hello,
>>
>> I'm working on implementing a simulation of a quantum computer in Sympy. One thing I have to do is make a matrix representation of a series of quantum gates and am having some troubles with the simplifications that automatically occur. I basically have two things that result in hard to read matrix expressions:
>>
>> 1. Sympy automatically turns exp(I*pi/2) into i and exp(I*pi) into -1, which results in ugly expressions like I*exp(I*pi/4) rather than the cleaner expression exp(3*I*pi/4).
>>
>> 2. I can't find a way to factor out constants that are common to all elements in a matrix. (In fact I can't even create a Mul object which holds a matrix and a constant because a sympy-matrix cannot be simpified)
>>
>> I've attached a screen-shot of one of the rather unpleasant looking matrices. Is there any way to fix either one of these problems?
>>
>> Thanks,
>> Addison J Cugini
>>
>
> --
> 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.
>
>
--
Brian E. Granger, Ph.D.
Assistant Professor of Physics
Cal Poly State University, San Luis Obispo
bgra...@calpoly.edu
elli...@gmail.com
Ronan Lamy
> Hello,
>
> I'm working on implementing a simulation of a quantum computer in
> Sympy. One thing I have to do is make a matrix representation of a
> series of quantum gates and am having some troubles with the
> simplifications that automatically occur. I basically have two things
> that result in hard to read matrix expressions:
>
> 1. Sympy automatically turns exp(I*pi/2) into i and exp(I*pi) into -1,
> which results in ugly expressions like I*exp(I*pi/4) rather than the
> cleaner expression exp(3*I*pi/4).
>
In [18]: exp(I*pi, evaluate=False)
Out[18]:
π⋅ⅈ
ℯ
However:
In [19]: 2*_
Out[19]: 2⋅(-1)
which looks like a bug to me.
> 2. I can't find a way to factor out constants that are common to all
> elements in a matrix. (In fact I can't even create a Mul object which
> holds a matrix and a constant because a sympy-matrix cannot be
> simpified)
These are deep and complex issues. Actually, I've been thinking about
matrices recently, so I have a few ideas.
First, I'll hijack your thread with my concerns: I've found two articles
([1], [2]) describing an algorithm for fraction-free LU decomposition.
The algorithm is already partly implemented as Matrix.LUdecompositionFF
but not really used. Briefly, given a matrix A with coefficients in some
integral domain I, the algorithm factors it as A = Pr^-1.L.D^-1.U.Pc
where Pr and Pc are permutation matrices, L is lower triangular, U is
upper triangular, D is diagonal and all matrices have their coefficients
in I. For symbolic or integer matrices, it happens to be asymptotically
faster than the conventional LU algorithm. Now, instead of storing the
result in 5 m*n matrices, it could be compacted into just one matrix
plus 2 permutations. This suggests creating some sort of LUDecomposition
class that has this compact internal representation, uses it for some
operations and knows how to expand back to the 5 full-fledged matrices
it represents.
Actually, there are many more things that could be done with specialised
classes. For instance, permutation matrices can be used much more
efficiently by using the underlying permutation: multiplying them with
another matrix is just permuting the rows or columns, their inverse is
just the matrix for the inverse permutation, etc. And there are many
types of matrices for which huge gains in efficiency are possible:
identity matrices, diagonal, orthogonal, hermitian,... Note that these
specialised matrix classes have to be immutable, since altering just one
element usually destroys the special property.
For all this to work, there needs to be an abstract base class for
matrices and, in particular, proper dispatch rules (see [3]) implemented
in __mul__/__rmul__ of each class.
Back to your problems. First, the only reasonable way to sympify a
Matrix is to make Matrix a subclass of Basic, which requires finding a
way to allow subclasses of Basic to be mutable.
Now, assume that you have an object like Mul(sqrt(2), 1/4, A), A being a
matrix. This object needs to know that it represents a matrix so that we
can multiply it by another matrix, take its inverse, etc. This would
imply that Mul is completely generic and defers to its computed type for
all attribute access, which doesn't really seem compatible with sympy's
current architecture.
So, sqrt(2)/4*A should be represented as a new class inheriting from the
matrix ABC. However, this is a heterogeneous product with a scalar and a
matrix so the structure should capture this asymmetry. A nice solution
could be something like LinearCombination({A: sqrt(2)/4}), which can
obviously also represent a*A+B as LinearCombination({A: a, B: 1}) and
behaves like just one more example of the specialised classes I
described above.
Finally, a somewhat unrelated remark: I find the Matrix API bloated and
quite inconsistent, I think it should be much closer to numpy and all
the methods duplicating functionality easily achieved with slicing
should be removed.
HTH,
Ronan
[1]:http://www.apmaths.uwo.ca/~djeffrey/Offprints/FFLUQR.pdf
[2]:http://www.apmaths.uwo.ca/~djeffrey/Offprints/David-Jeffrey-LU.pdf
[3]:http://docs.python.org/library/numbers.html#implementing-the-arithmetic-operations
Aaron S. Meurer
On Jul 13, 2010, at 4:08 PM, Ronan Lamy wrote:
> Le lundi 12 juillet 2010 à 11:41 -0700, Addison Cugini a écrit :
>> Hello,
>>
>> I'm working on implementing a simulation of a quantum computer in
>> Sympy. One thing I have to do is make a matrix representation of a
>> series of quantum gates and am having some troubles with the
>> simplifications that automatically occur. I basically have two things
>> that result in hard to read matrix expressions:
>>
>> 1. Sympy automatically turns exp(I*pi/2) into i and exp(I*pi) into -1,
>> which results in ugly expressions like I*exp(I*pi/4) rather than the
>> cleaner expression exp(3*I*pi/4).
>>
> You could do:
> In [18]: exp(I*pi, evaluate=False)
> Out[18]:
> π⋅ⅈ
> ℯ
>
> However:
>
> In [19]: 2*_
> Out[19]: 2⋅(-1)
>
> which looks like a bug to me.
evaluate=False expressions rarely work correctly after you do further things to them. I think maybe the whole thing needs refactoring (Mateusz's hold or something else).
I think it should be structured similar to the Polys. You have three levels. At the bottom level, it is just the pure Python representation (list or dictionary or whatever you use), and fast algorithms that work on just those, without any concern to container classes. Then, you have specialized classes to handle these, like the ones that you mentioned. These help separate implementations from the user (like sparse vs. dense, or special cases like permutations). These still require correct types in the arguments, and there isn't any automatic guessing. Then at the top, you have Matrix() (or Poly()), which is user friendly and makes decisions for the user if he doesn't specify. In Poly, this is things like what representation or domain to use; for Matrix, it would be things like what representation to use or if it should use one of your special, more efficient, classes.
We also need to work in the Poly ground types in the Matrices, as suggested in comment 6 of issue 1441.
>
> Back to your problems. First, the only reasonable way to sympify a
> Matrix is to make Matrix a subclass of Basic, which requires finding a
> way to allow subclasses of Basic to be mutable.
> Now, assume that you have an object like Mul(sqrt(2), 1/4, A), A being a
> matrix. This object needs to know that it represents a matrix so that we
> can multiply it by another matrix, take its inverse, etc. This would
> imply that Mul is completely generic and defers to its computed type for
> all attribute access, which doesn't really seem compatible with sympy's
> current architecture.
So we need to have issue 1941. By the way, it's not enough to just make __mul__ and __rmul__ work if you want to use it with symbols. Mul will always grab hold of certain things due to operator precedence and parentheses (trust me, I tried it last year with issue 1336).
> So, sqrt(2)/4*A should be represented as a new class inheriting from the
> matrix ABC. However, this is a heterogeneous product with a scalar and a
> matrix so the structure should capture this asymmetry. A nice solution
> could be something like LinearCombination({A: sqrt(2)/4}), which can
> obviously also represent a*A+B as LinearCombination({A: a, B: 1}) and
> behaves like just one more example of the specialised classes I
> described above.
>
> Finally, a somewhat unrelated remark: I find the Matrix API bloated and
> quite inconsistent, I think it should be much closer to numpy and all
> the methods duplicating functionality easily achieved with slicing
> should be removed.
I agree to make it close to numpy, which is what people will know. I don't see why we have to remove methods just because you can do it with a slice, though. Slicing can be confusing to new users, who may not even think to use it.
Aaron Meurer
smichr
On Jul 13, 12:26 am, "Aaron S. Meurer" <asmeu...@gmail.com> wrote:
> Well, fixing sympify to work with matrices would be a good start. There would probably be more work needed than that to make Matrix play nice with Mul, though.
>
> I think Chris had a branch somewhere with some kind of extract_gcd() method for Matrix. Maybe it could help you.
>
possible:
>>> A=Matrix([[sqrt(2)/4,-sqrt(2)/4]])
>>> A
[2**(1/2)/4, -2**(1/2)/4]
>>> A.gcdfactor()
2**(1/2)/4
>>> A
[1, -1]
>>> from sympy.abc import _EXP
>>> _EXP**(I*pi/2)*_EXP**(I*pi)
_E**(3*pi*I/2)
>>> _.subs(_EXP,S.Exp1)
-I
That branch is nearly error free and will soon be ready for review.
Brian Granger
I agree that this is not how most things are done in sympy. But as we
have been implementing the quantum stuff this summer, we are finding
that in some cases the current approach is quite difficult to work
with. For example, if we need to have custom Mul logic for our
operators, hilbert spaces or states, we have to add that logic to Mul
itself, which does not make sense. This is a separate issue from
being mutable, but I think in the long run, these issues will limit
the ability to sympy to grow and expand. Things to think about.
> So, sqrt(2)/4*A should be represented as a new class inheriting from the
> matrix ABC. However, this is a heterogeneous product with a scalar and a
> matrix so the structure should capture this asymmetry. A nice solution
> could be something like LinearCombination({A: sqrt(2)/4}), which can
> obviously also represent a*A+B as LinearCombination({A: a, B: 1}) and
> behaves like just one more example of the specialised classes I
> described above.
This is exactly the type of thing that we find ourselves wanting to
do. The problem is that it means that all, new non-trivial types will
end up implementing their own Mul class or subclass. This is smelly.
Cheers,
Brian
> Finally, a somewhat unrelated remark: I find the Matrix API bloated and
> quite inconsistent, I think it should be much closer to numpy and all
> the methods duplicating functionality easily achieved with slicing
> should be removed.
>
> HTH,
> Ronan
>
> [1]:http://www.apmaths.uwo.ca/~djeffrey/Offprints/FFLUQR.pdf
> [2]:http://www.apmaths.uwo.ca/~djeffrey/Offprints/David-Jeffrey-LU.pdf
> [3]:http://docs.python.org/library/numbers.html#implementing-the-arithmetic-operations
>
Brian Granger
>
>
> On Jul 13, 12:26 am, "Aaron S. Meurer" <asmeu...@gmail.com> wrote:
>> Well, fixing sympify to work with matrices would be a good start. There would probably be more work needed than that to make Matrix play nice with Mul, though.
>>
>> I think Chris had a branch somewhere with some kind of extract_gcd() method for Matrix. Maybe it could help you.
>>
>
> In the rev2 branch of smichr's github account the following is
> possible:
>
>>>> A=Matrix([[sqrt(2)/4,-sqrt(2)/4]])
>>>> A
> [2**(1/2)/4, -2**(1/2)/4]
>>>> A.gcdfactor()
> 2**(1/2)/4
>>>> A
> [1, -1]
Very nice, that would be quite useful for the quantum computing stuff.
>>>> from sympy.abc import _EXP
>>>> _EXP**(I*pi/2)*_EXP**(I*pi)
> _E**(3*pi*I/2)
>>>> _.subs(_EXP,S.Exp1)
> -I
>
> That branch is nearly error free and will soon be ready for review.
Great!
Brian
Ondrej Certik
What are the other options?
Ondrej
Brian Granger
The other option is strong encapsulation. By this, I mean that:
1. sympy types (subclasses of Basic) should contain methods and
attributes that implement *ALL* the logic related to that type.
2. Logic related to a sympy type should not be anywhere but its own
methods and classes.
With this type of encapsulation, sympy would not have logic that
switched on the type (if isinstance..., if is_*, etc.) that it
currently has throughout.
Applying this idea to Mul, we would have to go through Mul and
identify the types operations that Mul does (flatten, combining
powers, etc). Once we identify the generic things we need to do with
operands of Mul, we would move that logic out of Mul and into the
types themselves. Mul would simply call the appropriate methods of
its operands to accomplish type specific things. If types want to
work with Mul, they would need to implement the appropriate methods.
My intuition says that this type of design would work better with
mutable types.
Sympy follows this type of thing in certain places, but in other
places it doesn't. Doing this throughout sympy would be a massive
change though. But, as I create custom sympy types, this is the type
of thing that I continually run into (a lack of strong encapsulation).
I think for basic types it is less of an issue though.
Brian
> Ondrej
Aaron S. Meurer
The solution is the issue 1941 that I mentioned earlier. We need a Mul (and Add) that ask objects how to combine themselves, instead of forcing you to either modify Mul.flatten() or write your own incompatible Mul (like we worked on in Los Alamos last year). There are already several things in sympy right now that could be made simpler with this, such as Order and Infinity, but it would also help with things like units, non-commutative symbols, matrices, and other strange objects that can have scalar coefficients but that require special logic to "multiply". I think if we make a self-combination system that general and strong enough, there will be little need for separate Mul's.
It won't be easy to implement, but until it is, the core Mul will continue to be useless for anything that requires special multiplication rules without further special casing Mul.flatten(). (actually, there is a sample core implementation on the issue page, but merging in any change in the core is almost always a major headache).
By the way, this is what Pow already does with ._eval_power(), and it works very well. Of course, exponentiation is easier because it is non-commutative and non-associative.
Aaron Meurer
Brian Granger
Yes, this issue does address this exactly. I think it would be a very
nice thing to have.
> It won't be easy to implement, but until it is, the core Mul will continue to be useless for anything that requires special multiplication rules without further special casing Mul.flatten(). (actually, there is a sample core implementation on the issue page, but merging in any change in the core is almost always a major headache).
>
> By the way, this is what Pow already does with ._eval_power(), and it works very well. Of course, exponentiation is easier because it is non-commutative and non-associative.
Yes, it would be more difficult but the design pattern is the same. I
like how _eval_power is handled and I think we should do things like
that more.
Cheers,
Brian
> Aaron Meurer
Aaron S. Meurer
Exactly!
>
> With this type of encapsulation, sympy would not have logic that
> switched on the type (if isinstance..., if is_*, etc.) that it
> currently has throughout.
>
> Applying this idea to Mul, we would have to go through Mul and
> identify the types operations that Mul does (flatten, combining
> powers, etc).
Just a note, if it turns out that the standard combining of coefficients and exponents (x + x => 2*x, x*x => x**2) is slower this way (I don't know if it will be or not), we should leave that logic in Mul.flatten() for performance, because it is so fundamental for that particular operation.
> Once we identify the generic things we need to do with
> operands of Mul, we would move that logic out of Mul and into the
> types themselves. Mul would simply call the appropriate methods of
> its operands to accomplish type specific things. If types want to
> work with Mul, they would need to implement the appropriate methods.
> My intuition says that this type of design would work better with
> mutable types.
I don't see how they conflict (or even how mutable types would solve this particular problem).
Aaron Meurer
Brian Granger
The reason that I bring up mutable types is this. There are certain
types of operations that involve more than one type. For some of
these things, it may be difficult to isolate the logic into methods of
a single types. In situations like this, I think that it might make
more sense to just build the expression tree in place without doing
transformations, but then later going back and applying transformation
rules in place to the mutable expression tree.
Ondrej Certik
There is also the issue of speed. The basic manipulation has to be
fast. So we might have some speedy core, that is not so general and
then some slower abstraction, that allows you to do some cool tricks,
but it'll be slower.
In principle, I think we can have multiple things like that in sympy,
just like we have old and new assumptions. We'll see what works and
what doesn't and we can improve on the way.
Ondrej
Brian Granger
Yes, it may make sense to handle fundamental types in an optimized
manner. Definitely.
> In principle, I think we can have multiple things like that in sympy,
> just like we have old and new assumptions. We'll see what works and
> what doesn't and we can improve on the way.
Yep, I think that approach is serving sympy well.
Cheers,
Brian
> Ondrej
Andy Ray Terrel
I agree with Ronan here. I'm working on a application that turns
symbolic linear algebra routines into fast algorithms, right now I'm
rolling my own symbolic Matrices but would be interested in working on
a richer abstraction layer as well.
Essentially the current Matrix class is a container with some extra
routines. My preference would be to separate out the container view
from the algebraic structure view. Similar to numpy's array versus
matrix.
So is this yet another rewrite of the core?
The situation is still somewhat ugly, I would suggest the mul to have
some precedences to get a nice combination of things. For example,
when combining a constant (c) and piecewise function (p), does both
the constant.__mul__ and piecewise.__mul__ have to know about each
other? If they don't c*p != p*c . Perhaps with a precedence the
default mul will call the __mul__ of the other object then p*c -> p*c
(assuming c is higher than p).
asmeurer
Actually, I just noticed that sympy can reduce exp(3*I*pi/4) using
expand(complex=True):
In [1]: exp(3*I*pi/4)
Out[1]:
3⋅π⋅ⅈ
─────
4
ℯ
In [2]: exp(3*I*pi/4).expand(complex=True)
Out[2]:
⎽⎽⎽ ⎽⎽⎽
╲╱ 2 ⅈ⋅╲╱ 2
- ───── + ───────
2 2
So I think exp() should either always automatically reduce to radicals/
sin,cos for numeric arguments, or never (probably never is better).
Aaron Meurer