#32318 closed enhancement (fixed)

Generate non-redundant indices of tensors with symmetries efficiently

Reported by: Marius Gerbershagen Owned by:
Priority: major Milestone: sage-9.5
Component: performance Keywords:
Cc: Samuel Lelièvre, Michael Jung, Eric Gourgoulhon, Travis Scrimshaw Merged in:
Authors: Marius Gerbershagen Reviewers: Eric Gourgoulhon, Michael Jung
Report Upstream: N/A Work issues:
Branch: 08ccd3b (Commits, GitHub, GitLab) Commit: 08ccd3b89d25893230bbef55b0ed448b23c9d3de
Dependencies: Stopgaps:

Status badges

Description

In sage.tensor.modules.comp.CompWithSym, the method generating non-redundant indices is implemented very inefficiently by going through all indices and throwing away the redundant ones. In the worst case scenario of a rank n tensor in n dimensions, we have to iterate through nn indices only to find the single non-redundant one.

This ticket fixes this by changing the implementation to directly generate only non-redundant indices. This gives for example the following improvement for 7 dimensions:

sage: import timeit
sage: from sage.tensor.modules.comp import Components, CompWithSym, CompFullySym, CompFullyAntiSym
sage: V = VectorSpace(QQ, 7)
sage: c = CompFullyAntiSym(QQ, V.basis(), 7)
sage: timeit.timeit(lambda:list(c.non_redundant_index_generator()),number=1)

new implementation:
4.3510999603313394e-05

old implementation:
2.4149506480007403

I also simplified the method generating redundant indices a bit.

Change History (25)

comment:1 Changed 17 months ago by Marius Gerbershagen

Status: newneeds_review

comment:2 Changed 17 months ago by Marius Gerbershagen

Summary: Performace improvements in generating non-redundant indices of tensors with indicesPerformace improvements in generating non-redundant indices of tensors with symmetries

comment:3 Changed 17 months ago by Samuel Lelièvre

Cc: Samuel Lelièvre added
Summary: Performace improvements in generating non-redundant indices of tensors with symmetriesGenerate non-redundant indices of tensors with symmetries efficiently

comment:4 Changed 16 months ago by git

Commit: a60d93c37102245cb7f9e319576d9f92e7d15c1e1c108d85c05ccd171548f14f21b1c5a7db896d93

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

1c108d8tensor module components: non-redundant index generator: performance improvements

comment:5 Changed 16 months ago by Matthias Köppe

Cc: Michael Jung Eric Gourgoulhon added

comment:6 Changed 16 months ago by Michael Jung

Haven't checked the algorithm yet, but here are already some minor comments:

  • Why is self._sym copied but self._antisym not? From the first glimpse, a copy doesn't seem to be necessary anyway.
  • Very minor (PEP norm):
    - for k in range(1,len(isym)):
    + for k in range(1, len(isym)):
    

comment:7 Changed 16 months ago by Michael Jung

Very nice by the way. This could make a significant difference in sophisticated computations (especially since this generator is not cached). I am very curious about the timings of the SageManifolds? example notebooks!

Last edited 16 months ago by Michael Jung (previous) (diff)

comment:8 Changed 16 months ago by Eric Gourgoulhon

Status: needs_reviewneeds_work

Thanks for improving significantly the generator of non-redundant indices! There is an issue with the ordering of the output though, which affects the display of Christoffel symbols: the output shall respect the ordering of the coordinates, i.e. the doctests should not be modified.

comment:9 in reply to:  8 ; Changed 16 months ago by Marius Gerbershagen

Replying to egourgoulhon:

Thanks for improving significantly the generator of non-redundant indices! There is an issue with the ordering of the output though, which affects the display of Christoffel symbols: the output shall respect the ordering of the coordinates, i.e. the doctests should not be modified.

Is that really important? The order in which the components of tensors are displayed is a matter of personal taste, the mathematical results are the same anyway. Generating the indices in exactly the same order as we did previously would make the algorithm significantly more complicated.

comment:10 in reply to:  9 ; Changed 16 months ago by Eric Gourgoulhon

Replying to gh-spaghettisalat:

Replying to egourgoulhon:

Thanks for improving significantly the generator of non-redundant indices! There is an issue with the ordering of the output though, which affects the display of Christoffel symbols: the output shall respect the ordering of the coordinates, i.e. the doctests should not be modified.

Is that really important?

Yes I think so. Christoffel symbols are usually lengthy objects and displaying them in random order is not very useful. Ordered display is required for publishing or comparing with previously known results (e.g. textbook, wikipedia).

The order in which the components of tensors are displayed is a matter of personal taste, the mathematical results are the same anyway. Generating the indices in exactly the same order as we did previously would make the algorithm significantly more complicated.

Can one devise another method, just for display purposes, which would reorder the output of non_redundant_index_generator?

comment:11 in reply to:  10 Changed 16 months ago by Marius Gerbershagen

Replying to egourgoulhon:

Replying to gh-spaghettisalat:

Replying to egourgoulhon:

Thanks for improving significantly the generator of non-redundant indices! There is an issue with the ordering of the output though, which affects the display of Christoffel symbols: the output shall respect the ordering of the coordinates, i.e. the doctests should not be modified.

Is that really important?

Yes I think so. Christoffel symbols are usually lengthy objects and displaying them in random order is not very useful. Ordered display is required for publishing or comparing with previously known results (e.g. textbook, wikipedia).

The order is not random, it is just slightly different than before. For the christoffel symbols, previously the indices would increase from right to left, now the first index increases before the last two ones. In general, the new algorithm first increases unsymmetrized indices, then symmetrized and finally antisymmetrized ones. One could also do this in a different order (for instance first increase the antisymmetrized indices), which would restore the old order at least for christoffel symbols. But I don't really see the point of that, since there is no universally agreed standard for in which order one should display the indices.

The order in which the components of tensors are displayed is a matter of personal taste, the mathematical results are the same anyway. Generating the indices in exactly the same order as we did previously would make the algorithm significantly more complicated.

Can one devise another method, just for display purposes, which would reorder the output of non_redundant_index_generator?

You could of course just convert all of the indices into a list and sort that one via list(tensor.non_redundant_index_generator()).sort(). This is a bit inefficient but for display purposes, this should be fine.

comment:12 Changed 16 months ago by Michael Jung

Cc: Travis Scrimshaw added

I agree with Eric. Ordering the indices by coordinates is extremely convenient, at least for the output.

Since this seems only affect the Christoffel symbols, I can live with sorting them afterwards, even for the price of a small slow down in that case (which, I suspect, is still faster than the current implementation).

On the other hand, please keep in mind that there is a refactoring of the components's module going on in #30307. To that end, we will cache results related to symmetry such as this iterator.

Last edited 16 months ago by Michael Jung (previous) (diff)

comment:13 Changed 16 months ago by git

Commit: 1c108d85c05ccd171548f14f21b1c5a7db896d9308ccd3b89d25893230bbef55b0ed448b23c9d3de

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

08ccd3btensor module components: non-redundant index generator: performance improvements

comment:14 Changed 16 months ago by Marius Gerbershagen

I have implemented sorting of the indices for display. All outputs should now appear the same way as before.

Changed 6 hours ago by gh-mjungmath:

Haven't checked the algorithm yet, but here are already some minor comments:

  • Why is self._sym copied but self._antisym not? From the first glimpse, a copy doesn't seem to be necessary anyway.

Because this is modified later on. The new algorithm treats non-symmetrized indices as a symmetrization of a single index and adds those to the sym list.

  • Very minor (PEP norm):
        - for k in range(1,len(isym)):
        + for k in range(1, len(isym)):
    

This has been fixed.

Very nice by the way. This could make a significant difference in sophisticated computations (especially since this generator is not cached). I am very curious about the timings of the SageManifolds? example notebooks!

As long as these examples are in dimensions below 6, I doubt that this will make a significant impact.

comment:15 Changed 16 months ago by Marius Gerbershagen

Status: needs_workneeds_review

comment:16 Changed 16 months ago by Eric Gourgoulhon

Reviewers: Eric Gourgoulhon, Michael Jung
Status: needs_reviewpositive_review

Thanks for the changes! Michael, do you agree with the positive review?

comment:17 Changed 16 months ago by Michael Jung

I really dislike the usage of flags within loops, in this case step_finished. My intuition tells me, it should be possible to simplify the code much further. But I immediately don't see how.

But this is personal taste, and for now this should be fine.

comment:18 Changed 16 months ago by Michael Jung

This is a huge improvement though. Thank you Marius!

comment:19 in reply to:  17 ; Changed 16 months ago by Marius Gerbershagen

Replying to gh-mjungmath:

I really dislike the usage of flags within loops, in this case step_finished. My intuition tells me, it should be possible to simplify the code much further. But I immediately don't see how.

If only python had a goto statement...

comment:20 Changed 16 months ago by Samuel Lelièvre

Minor comments, feel free to ignore any or all of them.

Possible simplification for getting a sorted list:

-            generator = list(self.non_redundant_index_generator())
-            generator.sort()
+            generator = sorted(self.non_redundant_index_generator())

PEP8 recommends 2 spaces before inline comments (and some suggested rephrasings):

-                    return # end point reached
+                    return  # end point reached
-        sym = self._sym.copy() # we may modify this in the following
+        sym = self._sym.copy()  # make a copy we can modify
-            step_finished = False # each step generates a new index
+            step_finished = False  # each step generates a new index
+                    return # end point reach
+                    return  # end point reached

Longer comments can go on their own line(s) (and sometimes be made a bit shorter):

-                sym.append([pos]) # treat non-symmetrized indices as being symmetrized with themselves
+                # treat non-symmetrized indices as symmetrized with themselves
+                sym.append([pos])
-                    return # we went through all indices and didn't
-                           # find one which we can increase, thus we
-                           # have generated all indices
+                    # we went through all indices and none could be
+                    # increased, so we have generated all indices
+                    return

The following give the same boolean:

  • not any([a in b for b in c]) and not any([a in b for b in d])
  • all(a not in b for e in (c, d) for b in e)

so one could rewrite:

-            if not any([pos in isym for isym in sym]) and not any([pos in isym for isym in antisym]):
+            if all(pos not in isym for xsym in (sym, antisym) for isym in xsym):

Calling the flag done rather than step_finished, one would read: if not done: ....

Possible rewriting given how Python turns numbers and lists to booleans:

-                if not step_finished and i == 0 and len(antisym) == 0:
+                if not step_finished and not i and not antisym:

Maybe shorten this comment:

-                            # this index is at the maximum and we have
-                            # to reset it
+                            # this index is at the maximum; reset it

Maybe add spaces around + and -, and after commas, in several instances of:

- range(len(sym)-1,-1,-1)
+ range(len(sym) - 1, -1, -1)

comment:21 Changed 16 months ago by Michael Jung

I ponder, is the outer while-loop even necessary? Couldn't you just split this whole thing into two separated loops and yield within those? Then the code gets much simpler and the flag becomes obsolete.

Last edited 16 months ago by Michael Jung (previous) (diff)

comment:22 in reply to:  19 Changed 16 months ago by Michael Jung

Replying to gh-spaghettisalat:

Replying to gh-mjungmath:

I really dislike the usage of flags within loops, in this case step_finished. My intuition tells me, it should be possible to simplify the code much further. But I immediately don't see how.

If only python had a goto statement...

http://entrian.com/goto/ :P

comment:23 in reply to:  20 Changed 16 months ago by Eric Gourgoulhon

Replying to slelievre:

Possible rewriting given how Python turns numbers and lists to booleans:

-                if not step_finished and i == 0 and len(antisym) == 0:
+                if not step_finished and not i and not antisym:

Personally, I don't agree with that one. Of course, not i is more pythonic, but I find it less readable than i == 0, especially in the current context, where i is an integer.

comment:24 Changed 16 months ago by Matthias Köppe

Milestone: sage-9.4sage-9.5

comment:25 Changed 15 months ago by Volker Braun

Branch: public/tensor-indices08ccd3b89d25893230bbef55b0ed448b23c9d3de
Resolution: fixed
Status: positive_reviewclosed
Note: See TracTickets for help on using tickets.