tensors in SymPy

326 views
Skip to first unread message

mario

unread,
Jun 13, 2012, 7:38:55 AM6/13/12
to sy...@googlegroups.com
In https://github.com/sympy/sympy/pull/1326 it is implemented
the canonicalization of Riemann invariants using monoterm symmetries;
it is much faster than xperm.c from xPerm [1];
the reason is that in xperm.c there is duplication of tensors,
related to the intersection of the dummy symmetry group with a conjugate
of the slot symmetry group; this e.g. makes the difference between
O(N^3.4) and exponential behavior in a class of Riemann invariants,
N being the number of Riemann tensors.

From the point of view of representing tensors in SymPy, the above one
is an easy problem, since there is only one tensor type,
only one index type, algebraic operations are absent.

I think that now we should start to implement tensors in SymPy;
here is an outline of the program:

1) reorder products of tensors in a standard way,
so that they can be read as a single tensor T; this is explained
in a paper by Portugal [2]; this step is absent in Riemann invariants.

2) Obtains the symmetries of the tensor T from the symmetries
of the tensors composing it; for instance a tensor of the type
T = R R A A p q
where R is the riemann tensor, A is a symmetric tensor, p and q are
vectors, and the indices are contracted in some way, with some free indices;
tT inherits the slot symmetries of the single tensors, plus
the symmetries due to the fact that the R's commute; same for the S's.
This step is present in riemann.py with RiemannMonomial, and it is easily
generalized to more tensor types.
Then one canonicalizes T; this is done with canonicalize in
perm_algorithms.py

3) Add to T algebraic operations; canonicalization allows simplification.

4) Deal with multiterm symmetries, like the cyclic symmetry
of the Riemann tensor,

5) and with dimensional dependent identities, like the fact that a
sum of tensors which is antisymmetric in d+1 indices in d dimensions is zero.

For the last two steps probably one has to develop a module on Young tableaux;
I saw there is one in cadabra [3]

More points?

It seems to me that one can work fairly well on the various points separately.
I think that after implementing the first three points one can do already
a lot of interesting things.
Since canonicalization is very fast, I expect that tensors in SymPy can
be fast.

I would appreciate comments and I hope someone is interested in working on some
of these points.

[1] see  http://www.xact.es/index.html
[2] R. Portugal "An Algorithm to Simplify Tensor Expressions"
  Comput. Phys. Commun. 115, p.215 (1998)
[3] http://cadabra.phi-sci.com/

krastano...@gmail.com

unread,
Jun 13, 2012, 8:30:25 AM6/13/12
to sy...@googlegroups.com, Comer Duncan
=======

Just so we are sure that we speak about the same thing:

There are a few things that people mean when they say tensors (my
explanations and examples will include much hand-waving for the sake
of brevity):

1. multidimensional array supporting Einstein index notation (without
distinction between covariant and contravariant indices). This is
implemented (well I think) in the `Indexed` class in `tensors`.
However it is not of much use for anything other than code generation
for summation loops.


2. explicitly represented tensors in some basis. This is just a
multidimensional `Matrix` class. I have an ugly frankensteinian
implementation for such an object that I needed for a school project.
This is very close to the previous notion, however it supports
covariant/contravariant indices. Check
http://www.xact.es/xCoba/index.html for a very serious implementation
of this and more.


3. abstract tensors based on the index notation formalism (but with
Penrose abstract indices). This is what the core of xAct does, and
this is what Mario proposes (I think, correct me if I am wrong)


4. abstract tensors based on products of form fields and vector fields
(what my GSoC Differential Geometry project may provide (however it is
a secondary goal that may get finished after the end of the project))


An analogy to make things clearer:

1 relates to 2 relates to 3 relates to 4 as
numpy to matrix algebra to linear algebra to differential geometry

numpy - just arrays
matrix algebra - notions for certain high level operations
linear algebra - work independent of the choice of basis
differential geometry - vectors become vector fields, etc.

=======

I have also cc-ed this to prof. Comer Duncan as he has expressed
interest in using sympy for tensor-heavy workflow.

=======

My comments on the proposal:

One thing seems unclear to me: how will tensors be implemented. At the
moment this problem is apparent in other sympy modules:

- quantum mechanics - the need to represent kets and bras -
subclassing Expr and dealing with the need to implement the new logic
for handling the expression trees.

- matrix expression - the need to represent
abstract_matrix_symbol*matrix_that_actually_contains_elements, etc -
subclassing Expr and creating stuff like MatAdd, MatMul, etc and again
dealing with the need to implement new logic for handling the
expression trees.

A few times I have said that I do not think that this is sustainable.
What happens when I want to work with object from different
submodules... Instead I would prefer just to build the expression
trees using simple Add and Mul nodes and use crawlers that do what
they have to do. It would be much easier to add new logic to these
crawlers, like knowing what certain matrices can not be multiplied
with each other for instance. I am doing this for my Differential
Geometry project and this is also what the xAct project have done for
their tensor representation.

So I urge that when work starts on this project the author does not
create anything else than a simple Tensor class. He can then design
all the rest as functions that expect expressions trees containing
Tensor instances and other stuff in some canonical form. The
canonicalizer is separate. The function that verifies that the input
is not nonsense (Dividing scalar by a vector) is separate.

=======

One more thing:
> For the last two steps probably one has to develop a module on Young
> tableaux;
> I saw there is one in cadabra [3]
@Aleksandar Makelov, is this part of you GSoC project?

How much of this will depend on new functionality in the group theory module?

More generally, how can the work on canonicalization be done so it
does not interfere with the Aleksandar's work on group theory? It
would be a pity if you both implement the same stuff instead of
separating the work.

mario

unread,
Jun 13, 2012, 11:01:29 AM6/13/12
to sy...@googlegroups.com, Comer Duncan


On Wednesday, June 13, 2012 2:30:25 PM UTC+2, Stefan Krastanov wrote:
=======

....


3. abstract tensors based on the index notation formalism (but with
Penrose abstract indices). This is what the core of xAct does, and
this is what Mario proposes (I think, correct me if I am wrong)


4. abstract tensors based on products of form fields and vector fields
(what my GSoC Differential Geometry project may provide (however it is
a secondary goal that may get finished after the end of the project))


An analogy to make things clearer:

1 relates to 2 relates to 3 relates to 4 as
numpy to matrix algebra to linear algebra to differential geometry

numpy - just arrays
matrix algebra - notions for certain high level operations
linear algebra - work independent of the choice of basis
differential geometry - vectors become vector fields, etc.


abstract tensors also do not depend on the dimension; for instance the following
Riemann invariant (-d1 means covariant d1, d1 means contravariant d1)
is independent on the space-time dimension

```
>>> from riemann import RiemannMonomial
>>> s = 'R(-d4,-d2,-d6,-d3)*R(-d1,-d5,d2,d4)*R(d3,d6,d1,d5)'
>>> r = RiemannMonomial.from_string(s)
>>> r.canonic()
R(d1,d2,d3,d4)*R(-d1,-d2,d5,d6)*R(-d3,-d4,-d5,-d6)
```

 
One more thing:
> For the last two steps probably one has to develop a module on Young
> tableaux;
> I saw there is one in cadabra [3]
@Aleksandar Makelov, is this part of you GSoC project?

How much of this will depend on new functionality in the group theory module?

More generally, how can the work on canonicalization be done so it
does not interfere with the Aleksandar's work on group theory? It
would be a pity if you both implement the same stuff instead of
separating the work.

There is almost no interference; I have factored out the code for  orbit and orbit_traversal in
perm_groups.py, to be able to use them as functions without creating a PermutationGroup,
but that is all.
Aleksandar claims that his implementation of them is faster; I have not benchmarked them;
if they are it is great, and should be used everywhere (i.e. also in schreier_tree and in
the orbit_traversal function)

In the near future I do not plan to work on permutation groups, though I will follow with interest
Aleksandar's work, apart possibly from the following topic.

@Aleksandar Makelov do you plan to work on the intersection of groups?
That is somehow related to canonicalization, since the Butler algorithm for finding the
double coset representative is based on the intersection algorithm.

mario

unread,
Jun 13, 2012, 5:02:39 PM6/13/12
to sy...@googlegroups.com, Comer Duncan


On Wednesday, June 13, 2012 2:30:25 PM UTC+2, Stefan Krastanov wrote:
=======
....


 
I do not understand your remarks; maybe if I make a concrete example
you could tell me how you expect that algebraic operations should be implemented.
Suppose I write

```
t1^{a b c d} = g^{a b} p^c p^d
t2^{a b c d} = p^d p^c g^{b a}
```
where `p` is commuting and `g` is the metric.
I think that when the tensor objects are constructed they are put in canonical form;
do you agree with this?
Suppose that in the canonical form they are both given by the tensor object

```        heads, ordered indices, index positions,    slot symmetries
Tensor([g,p,p], [a,b,c,d],             [0,1,2,3],             [(0,1),(2,3)])
```
Now compute `t1^{a b c d} + t2^{a b c d}`;  which methods or functions
should I call to get  `2*g^{a b} p^c p^d` ?






Aleksandar Makelov

unread,
Jun 13, 2012, 6:21:22 PM6/13/12
to sy...@googlegroups.com
Yes, it seems that interference is not going to be a major problem.

> Aleksandar claims that his implementation of them is faster; I have not
> benchmarked them;
> if they are it is great, and should be used everywhere (i.e. also in
> schreier_tree and in
> the orbit_traversal function)

I tried both versions of orbit, stabilizer, orbit_transversal ( I
renamed it from orbit_traversal to orbit_transversal to keep with the
accepted terminology) on large symmetric/alternating/dihedral groups
and the backtrack searches were easily seen to perform worse. However,
there is some probability that they're more efficient on some other
types of groups, I don't know (I'll soon try to get some database of
known groups going, by the way). Could you provide some information on
their complexity, or some references about these algorithms? Or indeed
any literature related to computational group theory :) I've only been
using "Handbook of computational group theory" by Holt, and
"Permutation Groups" by Seress.

>
> In the near future I do not plan to work on permutation groups, though I
> will follow with interest
> Aleksandar's work, apart possibly from the following topic.
>
> @Aleksandar Makelov do you plan to work on the intersection of groups?
> That is somehow related to canonicalization, since the Butler algorithm
> for finding the
> double coset representative is based on the intersection algorithm.

There is an algorithm for finding the intersection of two subgroups of
a permutation group that uses some sort of backtracking, and it's
described in the "Handbook...", pp. 124-125. As far as I can see now,
Seress has only special cases of group intersections. Other than that,
I'm clueless. Do you have any resources on what algorithms for group
intersections might be efficient (I guess the one in the "Handbook..."
is not going to be terrible) ? I can go over it and implement it at
some point, how soon do you need it?
>
> --
> 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/-/BmRwhaM2irUJ.
>
> 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.

mario

unread,
Jun 14, 2012, 4:03:20 AM6/14/12
to sy...@googlegroups.com


On Thursday, June 14, 2012 12:21:22 AM UTC+2, Aleksandar Makelov wrote:



I tried both versions of orbit, stabilizer, orbit_transversal ( I
renamed it from orbit_traversal to orbit_transversal to keep with the
accepted terminology) on large symmetric/alternating/dihedral groups
and the backtrack searches were easily seen to perform worse. However,
there is some probability that they're more efficient on some other
types of groups, I don't know  (I'll soon try to get some database of
known groups going, by the way). Could you provide some information on
their complexity, or some references about these algorithms? Or indeed
any literature related to computational group theory :) I've only been
using "Handbook of computational group theory" by Holt, and
"Permutation Groups" by Seress.

 The orbit_traversal I wrote is depth-first, yours is breadth-first; to avoid problems
of incompatibility it would be better to use always the same.
I did not make a benchmark comparison, but I like yours better because it is shorter.
orbit_transversal seems a strange name to me; I believe that in search one
talks of  depth- and breadth-first traversals. About using a database, in testing I used a script to read from Atlas the generators, without saving them in a file; I attach it.

The Handbook by Holt is great; Seress seems more difficult; there is the book
by Butler; I found useful as an introduction
Kreher, Stinson "Combinatorial algorithms", which gives an overview to permutations,
partitions, graphs, but it is less advanced.
 
> In the near future I do not plan to work on permutation groups, though I
> will follow with interest
> Aleksandar's work, apart possibly from the following topic.
>
> @Aleksandar Makelov do you plan to work on the intersection of groups?
> That is somehow related to canonicalization, since the Butler algorithm
> for finding the
> double coset representative is based on the intersection algorithm.

There is an algorithm for finding the intersection of two subgroups of
a permutation group that uses some sort of backtracking, and it's
described in the "Handbook...", pp. 124-125. As far as I can see now,
Seress has only special cases of group intersections. Other than that,
I'm clueless. Do you have any resources on what algorithms for group
intersections might be efficient (I guess the one in the "Handbook..."
is not going to be terrible) ? I can go over it and implement it at
some point, how soon do you need it?

I do not know which implementation is efficient;
I read that it is non-polynomial, but that in practice it is fast.
It is great if you implement it; it would be useful to me now :)  but take your time.

 
t3.py

krastano...@gmail.com

unread,
Jun 14, 2012, 4:55:56 AM6/14/12
to sy...@googlegroups.com
> I do not understand your remarks; maybe if I make a concrete example
> you could tell me how you expect that algebraic operations should be
> implemented.
> Suppose I write
>
> ```
> t1^{a b c d} = g^{a b} p^c p^d
> t2^{a b c d} = p^d p^c g^{b a}
> ```
> where `p` is commuting and `g` is the metric.
> I think that when the tensor objects are constructed they are put in
> canonical form;
> do you agree with this?
> Suppose that in the canonical form they are both given by the tensor object
>
> ```        heads, ordered indices, index positions,    slot symmetries
> Tensor([g,p,p], [a,b,c,d],             [0,1,2,3],             [(0,1),(2,3)])
> ```
> Now compute `t1^{a b c d} + t2^{a b c d}`;  which methods or functions
> should I call to get  `2*g^{a b} p^c p^d` ?

Do you expect the user to work with an expression in non-canonical
form at any time? I mean: when the user has two canonically expressed
instances of a Tensor class and when he adds them together, do you
automagically get another canonically expressed Tensor instance or do
you get something like `Add(Tensor1, Tensor2)`?

I was expecting the second answer. If it is actually the second
answer, you will need to have a way to work with these expression
trees of Tensors and I tried to outline different possibilities,
mentioning the one that I like.

However I may have misunderstood your proposal. Maybe you do not
expect the creation of any user-accessible expression trees. I do not
know whether this is preferable or not at the moment?

I am not completely sure that I conveyed my thoughts clearly. Please,
do not hesitate to ask me for a better explanation.

mario

unread,
Jun 14, 2012, 6:53:30 AM6/14/12
to sy...@googlegroups.com


>Do you expect the user to work with an expression in non-canonical
form at any time? I mean: when the user has two canonically expressed
instances of a Tensor class and when he adds them together, do you
automagically get another canonically expressed Tensor instance or do
you get something like `Add(Tensor1, Tensor2)`?

I was expecting the second answer. If it is actually the second
answer, you will need to have a way to work with these expression
trees of Tensors and I tried to outline different possibilities,
mentioning the one that I like.

 
The canonicalization is on tensor monomials; when one of them is created
it should be put in canonical form.
Two tensors are equal if they have the same canonical form.
So in general if  `T1` and `T2` are tensor monomials, their sum is
`Add(T1, T2)` unless they are equal, in which case it is `2*T1`

In seems to me that you think that `TensorMonomial` should not inherit from Expr or QExpr;


> Instead I would prefer just to build the expression
trees using simple Add and Mul nodes and use crawlers that do what
they have to do.

How `Add(T1, T2)` should be implemented then?


 

krastano...@gmail.com

unread,
Jun 14, 2012, 9:02:41 AM6/14/12
to sy...@googlegroups.com
I suggested that instead of a class akin to TensorMonomial you simply
use Add. But I do not know the details, so I may be wrong.

Aleksandar Makelov

unread,
Jun 15, 2012, 8:42:12 AM6/15/12
to sy...@googlegroups.com


14 юни 2012, четвъртък, 11:03:20 UTC+3, mario написа:


On Thursday, June 14, 2012 12:21:22 AM UTC+2, Aleksandar Makelov wrote:



I tried both versions of orbit, stabilizer, orbit_transversal ( I
renamed it from orbit_traversal to orbit_transversal to keep with the
accepted terminology) on large symmetric/alternating/dihedral groups
and the backtrack searches were easily seen to perform worse. However,
there is some probability that they're more efficient on some other
types of groups, I don't know  (I'll soon try to get some database of
known groups going, by the way). Could you provide some information on
their complexity, or some references about these algorithms? Or indeed
any literature related to computational group theory :) I've only been
using "Handbook of computational group theory" by Holt, and
"Permutation Groups" by Seress.

 The orbit_traversal I wrote is depth-first, yours is breadth-first; to avoid problems
of incompatibility it would be better to use always the same.
I did not make a benchmark comparison, but I like yours better because it is shorter.
orbit_transversal seems a strange name to me; I believe that in search one
talks of  depth- and breadth-first traversals. About using a database, in testing I used a script to read from Atlas the generators, without saving them in a file; I attach it.


Thanks a lot for the script! I'll compare the running times for both versions of orbit, orbit_transversal, and stabilizer (I think these are all the ones I've changed). What do you use for benchmarking? (I use ipython's %timeit but maybe that's sort of an amateur technique).
As for the naming convention, transversal is the name used in the "Handbook" and GAP and by Seress. I think it derives from the following definition:

In combinatorial mathematics, given a collection C of sets, a transversal is a set containing exactly one element from each member of the collection.

( http://en.wikipedia.org/wiki/Transversal_%28combinatorics%29 )

in the sense that, in general, there is more than one element that sends alpha to some point in the orbit of alpha.


 
The Handbook by Holt is great; Seress seems more difficult; there is the book
by Butler; I found useful as an introduction
Kreher, Stinson "Combinatorial algorithms", which gives an overview to permutations,
partitions, graphs, but it is less advanced.

Yep, I agree with what you say about the Handbook and the book by Seress :) I found "Combinatorial Algorithms", and I'm still trying to find the
one by Butler, but that's not a big issue right now.
 
> In the near future I do not plan to work on permutation groups, though I
> will follow with interest
> Aleksandar's work, apart possibly from the following topic.
>
> @Aleksandar Makelov do you plan to work on the intersection of groups?
> That is somehow related to canonicalization, since the Butler algorithm
> for finding the
> double coset representative is based on the intersection algorithm.

There is an algorithm for finding the intersection of two subgroups of
a permutation group that uses some sort of backtracking, and it's
described in the "Handbook...", pp. 124-125. As far as I can see now,
Seress has only special cases of group intersections. Other than that,
I'm clueless. Do you have any resources on what algorithms for group
intersections might be efficient (I guess the one in the "Handbook..."
is not going to be terrible) ? I can go over it and implement it at
some point, how soon do you need it?

I do not know which implementation is efficient;
I read that it is non-polynomial, but that in practice it is fast.
It is great if you implement it; it would be useful to me now :)  but take your time.

Ok, I'll try to implement the group intersection algorithm described in the "Handbook" in the next few days :)

 
Reply all
Reply to author
Forward
0 new messages