Location for Tr operation

[Guru Devanla]

Hello All,

Recently I have been working on the Trace operator and related operations including partial traces. Currently, this operator is used exclusively by objects in the quantum module, though the operation itself is generic.  Currently, this module resides in /sympy/core/trace.py  There were some review comments by @asmeurer and @flacjacket proposing different locations for this module. Some comments can be found at https://github.com/sympy/sympy/pull/1396#issuecomment-6715963.

An alternative place for this module could be  : sympy/matrices/expressions. Currently, this folder has other matrix related expressions such as  matadd, matmul, matpow, inverse and transpose.

Where do you suggest this should go?

-Guru

[Aaron Meurer]

That sounds like the best place for it to go, assuming it plays well with MatrixExpr objects. 

Aaron Meurer

Where do you suggest this should go?

-Guru

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To view this discussion on the web visit https://groups.google.com/d/msg/sympy/-/Rn94--N-zG8J.
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.

[Guru Devanla]

Aaron,

Right now Tr class inherits from Expr directly. I will have to look into MatExpr as base class if that is necessary.

-Guru

To unsubscribe from this group, send email to sympy+unsubscribe@googlegroups.com.

[Matthew Rocklin]

I think that it makes sense to include a generic Trace object in matrix expressions. I think that it should inherit from Expr and not MatExpr (in my mind trace is a function :: Matrix -> Scalar). 

Guru, if you want a code example I would emulate the current Transpose object. 

I think that the Trace object in matrices/expressions should be very minimal however and should rely on _eval_trace heavily. sympy.physics and sympy.matrices handle things very differently. If we're going to share an object then that object should be very general (I think it should assume less even than the current implementation of Transpose which I will fix in a few minutes). 

To view this discussion on the web visit https://groups.google.com/d/msg/sympy/-/PgabKT0NVaIJ.

To post to this group, send email to sy...@googlegroups.com.

To unsubscribe from this group, send email to sympy+un...@googlegroups.com.

[Aaron Meurer]

On Jul 12, 2012, at 10:26 AM, Matthew Rocklin <mroc...@gmail.com> wrote:

I think that it makes sense to include a generic Trace object in matrix expressions. I think that it should inherit from Expr and not MatExpr (in my mind trace is a function :: Matrix -> Scalar). 

Guru, if you want a code example I would emulate the current Transpose object. 

I think that the Trace object in matrices/expressions should be very minimal however and should rely on _eval_trace heavily. sympy.physics and sympy.matrices handle things very differently. If we're going to share an object then that object should be very general (I think it should assume less even than the current implementation of Transpose which I will fix in a few minutes). 

So do you think that the same object should trace over MatrixExpr, Matrix, and the physics classes?  It seems to me that it should. 

Aaron Meurer

[Matthew Rocklin]

So do you think that the same object should trace over MatrixExpr, Matrix, and the physics classes?  It seems to me that it should. 

Yes. I think that we can easily write a sufficiently general version of Trace, something like the Transpose object in this branch

https://github.com/mrocklin/sympy/blob/transpose-decentralize/sympy/matrices/expressions/transpose.py

[Brian Granger]

One thing to be aware of is that the general Trace has to be able to
handle partial traces of tensor products. This logic is implemented
in this branch:

https://github.com/sympy/sympy/pull/1396

And is somewhat specialized to the stuff in quantum, although it
doesn't *depend* on quantum.

Cheers,

Brian

--
Brian E. Granger
Cal Poly State University, San Luis Obispo
bgra...@calpoly.edu and elli...@gmail.com

[Matthew Rocklin]

I had to look up partial trace (this is a new concept for me). After reading this I'm inclined to separate Trace and PartialTrace classes in an effort to optimize/simplify for the common case (at least the common case as I see it). 

The idea of a trace of a tensor product still sounds strange to me. I'm inclined to call this a tensor contraction. As I said though I'm not familiar with this specific notion of trace.

[Brian Granger]

On Thu, Jul 12, 2012 at 2:19 PM, Matthew Rocklin <mroc...@gmail.com> wrote:
> I had to look up partial trace (this is a new concept for me). After reading
> this I'm inclined to separate Trace and PartialTrace classes in an effort to
> optimize/simplify for the common case (at least the common case as I see
> it).

Yes, the concept is quite subtle and not used (as far as I know)
outside the context of quantum mechanics (where it is used a lot). I
thought about having separate Trace and PartialTrace classes, but from
the user API standpoint, they are really the same thing.

Trace = Partial Trace over all indices

Separating them adds complexity and confusion. Also in the current
implementation, the partial trace logic doesn't affect the regular
trace logic.

> The idea of a trace of a tensor product still sounds strange to me. I'm
> inclined to call this a tensor contraction. As I said though I'm not
> familiar with this specific notion of trace.

It is sort of like a tensor contraction, but I don't think they are
the same thing.

[Matthew Rocklin]

Yes, the concept is quite subtle and not used (as far as I know)
outside the context of quantum mechanics (where it is used a lot).  I
thought about having separate Trace and PartialTrace classes, but from
the user API standpoint, they are really the same thing.

Trace = Partial Trace over all indices

The immediate issue I see if this goes into MatrixExprs is that PartialTrace is a MatExpr while Trace is an Expr/Scalar. It's not easy to have a class that is both at the same time. If we do this then my approach would be to make both classes and have PartialTrace(arg, indices) yield a new Trace object on object creation when indices == all_indices. 

However, if partial traces are rarely used outside of quantum then maybe it makes sense to just keep them in quantum. A lot of the code in Guru's PR doesn't seem to apply to the MatrixExpr case (namely the tensor product stuff). MatrixExpressions is pretty restricted to linear algebra and not multi-linear algebra. I'm inclined to keep it this way. I think that an attempt to generalize the current module will result in an ugly solution.

I think that this stuff should probably live in quantum until we get around to a Tensor expressions module.

I think that I should make a simple Trace for MatrixExpressions and that the more sophisticated partial-trace on tensor-products code should move from sympy/core to sympy/physics/quantum for now.

[Guru Devanla]

Based on the discussion so far, it seems that the opinion is current Tr module should be moved to sympy/physics/quantum. Matt, per you suggestion do you plan to add a simple Trace under MatrixExpressions, in that case?

Also, can I leave the name of the module the way it is now or should we change to name to reflect that it belongs to the physics/quantum module.

I plan to work on this item this week and wrap it up based on whatever decision is made.

Regards

Guru

[Aaron Meurer]

On Tue, Jul 31, 2012 at 4:55 PM, Guru Devanla <grd...@gmail.com> wrote:
> Based on the discussion so far, it seems that the opinion is current Tr
> module should be moved to sympy/physics/quantum. Matt, per you suggestion do
> you plan to add a simple Trace under MatrixExpressions, in that case?
>
> Also, can I leave the name of the module the way it is now or should we
> change to name to reflect that it belongs to the physics/quantum module.

If you're worried about compatibility, changing the path breaks just
as much as changing the file name, so you might as well change it
(both will break exact path imports).

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To view this discussion on the web visit

> https://groups.google.com/d/msg/sympy/-/APU5J_H7v3UJ.

[Guru Devanla]

I went through the entire thread again, I feel we should just move the 

1.  Tr module into matrices/expressions

2.  Decide between inheriting from Expr or MatExpr. I would let it inherit from MatExpr, but set is_Matrix = false, if indices = all_indices.

3.  And *not* have 2 classes, since except for indices all other implementation details are same, and user gets cleaner API to deal with.

Here is the summary of some points made and why I feel this way:

1. Tr == Partial trace on all indices.  This just says the arguments of Tr

   behave differently during evaluation based on the implementation of

   _eval_args(). Even if we had 2 different classes(Tr, PartialTr), the

   implementations would be the same since we end up calling

   _eval_trace(). Every other implementation details of these 2 classes would

   be the same.

2. The current implementation though is used heavily by quantum module, still

   does not *depend* on quantum module.  This just states that  Tr

   is not *required* to reside in the physics/quantum folder.

3. Matthew suggested we could have 2 different classes and let Partial Trace

   yield a Tr object if indices == all_indices.  But, looking from user API

   perspective this just complicates things. I was looking at the use cases

   when user would be creating these objects, and did not think we should

   force user to pick between Tr and PartialTr. It should just be Tr(args,

   [tuple of indices]).

   The one point Matthew made was regarding during evaluation Tr yields an

   Expr/Scalar whereas partial trace yields another matrix. But, is this point

   alone a good reason to split these classes, given all other implementation

   details are the same.

Alternatively, if we decide to split these classes and place both of them in

matrices/expressions. Then, I see that  the PartialTr class does not have

any specific functionality except to yield a Tr object in some cases. Apart

from that, all other implementation details on Trace class applies to PartialTr.

So, in summary I would move the Tr class to matrices/expressions and let in

inherit from MatExpr. The user could always use Tr() class to run either trace or

partial trace operation. During _eval_trace(), if the operation is a partial

trace operation, then an appropriate object ( matrix ) is returned.

Also, if indices = all_indices, we could set the is_Matrix = False.

Or, if MatExpr does not apply to general case of Tr, then we should simple let

Tr inherit from Expr. I am not sure if this inheritance heirarchy should stop

Tr from being in the matrices/expressions folder.

Thoughts?

-Guru

On Thursday, July 12, 2012 5:34:22 PM UTC-5, Matthew wrote:

[Brian Granger]

Guru,

Thanks for summarizing things and making a recommendation. I agree
with your proposal.

On Wed, Aug 1, 2012 at 9:33 AM, Guru Devanla <grd...@gmail.com> wrote:
> I went through the entire thread again, I feel we should just move the
>
> 1. Tr module into matrices/expressions
> 2. Decide between inheriting from Expr or MatExpr. I would let it inherit
> from MatExpr, but set is_Matrix = false, if indices = all_indices.
> 3. And *not* have 2 classes, since except for indices all other
> implementation details are same, and user gets cleaner API to deal with.
>
> Here is the summary of some points made and why I feel this way:
>
> 1. Tr == Partial trace on all indices. This just says the arguments of Tr
> behave differently during evaluation based on the implementation of
> _eval_args(). Even if we had 2 different classes(Tr, PartialTr), the
> implementations would be the same since we end up calling
> _eval_trace(). Every other implementation details of these 2 classes
> would
> be the same.

Yes, because of this I strongly feel Tr and PartialTr should not be
separate classes.

> 2. The current implementation though is used heavily by quantum module,
> still
> does not *depend* on quantum module. This just states that Tr
> is not *required* to reside in the physics/quantum folder.

Yep.

> 3. Matthew suggested we could have 2 different classes and let Partial Trace
> yield a Tr object if indices == all_indices. But, looking from user API
> perspective this just complicates things. I was looking at the use cases
> when user would be creating these objects, and did not think we should
> force user to pick between Tr and PartialTr. It should just be Tr(args,
> [tuple of indices]).

Again, another point for these being the same class.

> The one point Matthew made was regarding during evaluation Tr yields an
> Expr/Scalar whereas partial trace yields another matrix. But, is this
> point
> alone a good reason to split these classes, given all other
> implementation
> details are the same.

I don't have a problem with this. If this becomes an issue in general
matrix expressions, we can use the is_Matrix attribute to allow
determination of its return value.

>
> Alternatively, if we decide to split these classes and place both of them in
> matrices/expressions. Then, I see that the PartialTr class does not have
> any specific functionality except to yield a Tr object in some cases. Apart
> from that, all other implementation details on Trace class applies to
> PartialTr.
>
>
> So, in summary I would move the Tr class to matrices/expressions and let in
> inherit from MatExpr. The user could always use Tr() class to run either
> trace or
> partial trace operation. During _eval_trace(), if the operation is a partial
> trace operation, then an appropriate object ( matrix ) is returned.
>
> Also, if indices = all_indices, we could set the is_Matrix = False.

Yep.

> Or, if MatExpr does not apply to general case of Tr, then we should simple
> let
> Tr inherit from Expr. I am not sure if this inheritance heirarchy should
> stop
> Tr from being in the matrices/expressions folder.
>
>
> Thoughts?

I am +1 on this.

Brian

> -Guru
>
>
>
>
>
>
> On Thursday, July 12, 2012 5:34:22 PM UTC-5, Matthew wrote:
>>>
>>> Yes, the concept is quite subtle and not used (as far as I know)
>>> outside the context of quantum mechanics (where it is used a lot). I
>>> thought about having separate Trace and PartialTrace classes, but from
>>> the user API standpoint, they are really the same thing.
>>>
>>> Trace = Partial Trace over all indices
>>
>>
>> The immediate issue I see if this goes into MatrixExprs is that
>> PartialTrace is a MatExpr while Trace is an Expr/Scalar. It's not easy to
>> have a class that is both at the same time. If we do this then my approach
>> would be to make both classes and have PartialTrace(arg, indices) yield a
>> new Trace object on object creation when indices == all_indices.
>>
>> However, if partial traces are rarely used outside of quantum then maybe
>> it makes sense to just keep them in quantum. A lot of the code in Guru's PR
>> doesn't seem to apply to the MatrixExpr case (namely the tensor product
>> stuff). MatrixExpressions is pretty restricted to linear algebra and not
>> multi-linear algebra. I'm inclined to keep it this way. I think that an
>> attempt to generalize the current module will result in an ugly solution.
>>
>> I think that this stuff should probably live in quantum until we get
>> around to a Tensor expressions module.
>>
>> I think that I should make a simple Trace for MatrixExpressions and that
>> the more sophisticated partial-trace on tensor-products code should move
>> from sympy/core to sympy/physics/quantum for now.

[Matthew Rocklin]

I am still -1 on this. 

I have no problem to an unevaluated PartialTrace object living in matrix/expressions. I just think that it should be separate.

Here are some of my concerns. 

From a user perspective having two separate classes is simpler in expectation. 99% of the time a user will want Trace. These users should not have to deal with the complexity of PartialTrace. I.e. the docstring for Trace should be very simple and intuitive. PartialTrace is foreign to most users. 

Trace is a scalar, PartialTrace can be a Matrix. The type system in the MatrixExpr module is much stricter than in physics/quantum. There is no way to have an object that is both. All MatrixExpr's have a shape. What is the shape of Trace(X)? If it is 1x1 then it will not be able to interact with other MatrixExprs. The MatrixExpr module is very clean and simple right now. I will strongly object to any type magic. The module is currently very simple, I would really like to keep it this way.

I think that having two simple classes is far simpler than having one complex class that can be many things.

Here is a PR for simple Trace. It is on top of another PR so only the last two commits are relevant to this conversation

https://github.com/sympy/sympy/pull/1456

[Matthew Rocklin]

What about having a Trace(Expr) object, a PartialTrace(MatrixExpr) object (which must be strictly a partial trace, not a full one, I.e. indices can not equal all_indices) and a function that provides a general interface. 

I would be much happier with this. I would still push for 

1) _eval_trace belongs to Trace, _eval_partial_trace would be created for whatever objects want to use it (stuff in quantum)

2) the interface function is not the default way in which people interact with Trace. It is an optional convenience function for those who want to sometimes use PartialTrace. 

In this setup it is obvious what type each class is. There is a unified interface for those who want to use indices=....

[Brian Granger]

Matthew,

On Wed, Aug 1, 2012 at 10:47 AM, Matthew Rocklin <mroc...@gmail.com> wrote:
> I am still -1 on this.

There is a larger discussion lurking here though. Right now, the
vision of matrix expressions leaves out a huge class of operation
involving matrices, namely those related to tensor products. Many
operations on matrices naturally extend to tensor products of matrices
and all of these things also extend to linear operators on Hilbert
spaces. Other examples of things that can be defined on tensor
products of matrices: Transpose, Adjoint. All of these operations
have a similar API as Trace/Partial trace in that you can do the
operation on all matrices in a tensor product, or just some of them
(the "partial" version). I agree these operations are not common
outside of quantum mechanics, but tensor products are a mainstream
part of linear algebra.

Are you open to the matrix expression modules gaining knowledge of
tensor products? If not, we are going to have to reimplement
everything in both matrix expressions and outside of it. BTW, we are
already going in this direction, in the sense that we have a PR open
#1158 for transpose/adjoint classes in
sympy/functions/elementary/complexes.py. The quantum module is using
those. Maybe this is the direction we need to go. In that case, why
not put Trace in the same place as transpose/adjoint, as they are very
similar.

> I have no problem to an unevaluated PartialTrace object living in
> matrix/expressions. I just think that it should be separate.
>
> Here are some of my concerns.
>
> From a user perspective having two separate classes is simpler in
> expectation. 99% of the time a user will want Trace. These users should not
> have to deal with the complexity of PartialTrace. I.e. the docstring for
> Trace should be very simple and intuitive. PartialTrace is foreign to most
> users.

Unless a user is doing quantum mechanics, then they will want partial
traces a good fraction of the time. I am still thinking about having
separate trace and partial trace operations, but I want to come up
with a unified way of handling these partial operations (for
transpose/adjoint) as well. Not sure it makes sense to have partial
transpose/adjoint as well. I am still leaning towards having a single
object, but am considering the 2 object option as well.

> Trace is a scalar, PartialTrace can be a Matrix. The type system in the
> MatrixExpr module is much stricter than in physics/quantum. There is no way
> to have an object that is both. All MatrixExpr's have a shape. What is the
> shape of Trace(X)? If it is 1x1 then it will not be able to interact with
> other MatrixExprs. The MatrixExpr module is very clean and simple right now.
> I will strongly object to any type magic. The module is currently very
> simple, I would really like to keep it this way.

I can see how this would be useful in matrix expressions if you don't
ever plan on generalizing things to include tensor products. If you
do though, these things will come up.

> I think that having two simple classes is far simpler than having one
> complex class that can be many things.

But it is not really that complex and it is not "many" things, just two at most.

> Here is a PR for simple Trace. It is on top of another PR so only the last
> two commits are relevant to this conversation
> https://github.com/sympy/sympy/pull/1456

I guess the larger question is how you see matrix expressions
co-existing with the rest of sympy. Is the stuff in matrix
expressions so specialized that it is justified to have its own
implementations of everything like trace/transpose/adjoint?

Cheers,

Brian

> --
> 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 Granger]

On Wed, Aug 1, 2012 at 11:02 AM, Matthew Rocklin <mroc...@gmail.com> wrote:
> What about having a Trace(Expr) object, a PartialTrace(MatrixExpr) object
> (which must be strictly a partial trace, not a full one, I.e. indices can
> not equal all_indices) and a function that provides a general interface.

That is a possibility, but I first want to understand the bigger
picture of where things like this are going to live in sympy and if
there will be duplication.

> --
> 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.

[Matthew Rocklin]

I'll reply more substantially in a bit (busy at work now) but I'll put in a quick word now. 

I think that an optimal version of SymPy would have a dedicated tensor module. I think that MatrixExpr would subclass TensorExpr. 

Matrix expressions would inherit a lot from tensor expressions but I still think that they would need to be part of a separate submodule. Linear algebra has been developed much more than multilinear algebra and there are some thoughts and many use-cases that don't generalize. 

I don't think that MatrixExprs should be generalized to TensorExprs. I think that TensorExprs should be created and then code from MatrixExprs factored out. 

(I haven't mentioned anything here about the Tensor vs NDArray discussion, which is also relevant but complicates this idea)

I'll take a look at your PR. Thanks for bringing it to my attention :)

[Brian Granger]

On Wed, Aug 1, 2012 at 1:55 PM, Matthew Rocklin <mroc...@gmail.com> wrote:
> I'll reply more substantially in a bit (busy at work now) but I'll put in a
> quick word now.
>
> I think that an optimal version of SymPy would have a dedicated tensor
> module. I think that MatrixExpr would subclass TensorExpr.

I should clarify something. The tensor product of matrices is itself
a matrix, not a higher dimensional tensor. By TensorExpr, I think you
are referring to something more general = a full blown tensor with
rank potentially > 2. Those types of tensors are not in any way
relevant in quantum mechanics.

> Matrix expressions would inherit a lot from tensor expressions but I still
> think that they would need to be part of a separate submodule. Linear
> algebra has been developed much more than multilinear algebra and there are
> some thoughts and many use-cases that don't generalize.

Again, multilinear algebra is not relevant here, just plain old linear
algebra.

> I don't think that MatrixExprs should be generalized to TensorExprs. I think
> that TensorExprs should be created and then code from MatrixExprs factored
> out.

I think this vision of having TensorExpr and MatrixExpr does make
sense, but again, the types of operations under discussion here are
under the matrix/linear algebra category.

[Matthew Rocklin]

Can you define tensor product then? Are you referring to hadamard / elementwise product?

[Brian Granger]

For matrices it looks like this:

http://en.wikipedia.org/wiki/Kronecker_product

Cheers,

Brian

[Matthew Rocklin]

Ah, I see, I would suggest that we call this the Kronecker product, not Tensor product. In my mind tensor product means this

http://en.wikipedia.org/wiki/Tensor_product

and more specifically this

http://en.wikipedia.org/wiki/Tensor_product#Tensor_product_of_two_tensors

but the term tensor is unfortunately overloaded. 

To me, Kronecker product, partial trace, etc... are all very obscure functionality. This is probably because I think about numerical linear algebra rather than quantum mechanics and because I am unfamiliar with them. My perspective is skewed. 

I am not in general opposed to anything joining MatrixExprs. I am opposed to these things causing what I think of as the MatrixExpr core to become more complex. I believe that this tension is relieved by having many atomic classes rather than a few swiss-army-knife classes. I ascribe fairly strongly to the idea that two simple classes are better than one complex class. 

Regarding the broader issue sympy/physics/quantum works very differently from MatrixExpr - at least that was the impression I gained when I studied it while building MatrixExprs. Quantum builds itself off of the SymPy core objects (i.e. Symbol, Expr) while MatrixExpr created new core objects (MatrixSymbol, MatrixExpr). I suspect that an attempt to glom together the two systems will result in tension. 

I care strongly about the core of MatrixExprs and have thoughts about exactly what it should look like and how it should grow. I do not however own MatrixExprs and, as far as I can tell, I'm currently the only user. So if it should evolve away from my vision of it (which looks like this) then that's ok. I'll probably work on it much less, but that's ok too.

If you're actually thinking of re-implementing everything in MatrixExprs then my proposed solution is to have a MatrixExpr core and have separate MatrixExpr modules. I would put PartialTrace into a module and once it proves useful for other MatrixExpr tasks we merge it into the core. If we wanted to merge Quantum and MatrixExprs I think this would also be a good approach. This is likely a lot of work though. 

It is however, I think, a good idea in the long run. There is work that I have planned for MatrixExprs (like an assumptions system) that I suspect could benefit quantum. The idea of PR1158 was to put a hermitian assumption into the SymPy core assumptions system. I think this sort of addition belongs in a new assumptions system for MatrixExprs. One that looks and operates much like this.

[Matthew Rocklin]

I'd like to merge in this PR for a Trace object in MatrixExprs. Any objections?

https://github.com/sympy/sympy/pull/1456