Opened 17 months ago
Closed 15 months ago
#32318 closed enhancement (fixed)
Generate nonredundant indices of tensors with symmetries efficiently
Reported by:  Marius Gerbershagen  Owned by:  

Priority:  major  Milestone:  sage9.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: 
Description
In sage.tensor.modules.comp.CompWithSym
, the method generating nonredundant 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 n^{n} indices only to find the single nonredundant one.
This ticket fixes this by changing the implementation to directly generate only nonredundant 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.3510999603313394e05 old implementation: 2.4149506480007403
I also simplified the method generating redundant indices a bit.
Change History (25)
comment:1 Changed 17 months ago by
Status:  new → needs_review 

comment:2 Changed 17 months ago by
Summary:  Performace improvements in generating nonredundant indices of tensors with indices → Performace improvements in generating nonredundant indices of tensors with symmetries 

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

Summary:  Performace improvements in generating nonredundant indices of tensors with symmetries → Generate nonredundant indices of tensors with symmetries efficiently 
comment:4 Changed 16 months ago by
Commit:  a60d93c37102245cb7f9e319576d9f92e7d15c1e → 1c108d85c05ccd171548f14f21b1c5a7db896d93 

comment:5 Changed 16 months ago by
Cc:  Michael Jung Eric Gourgoulhon added 

comment:6 Changed 16 months ago by
Haven't checked the algorithm yet, but here are already some minor comments:
 Why is
self._sym
copied butself._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
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!
comment:8 followup: 9 Changed 16 months ago by
Status:  needs_review → needs_work 

Thanks for improving significantly the generator of nonredundant 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 followup: 10 Changed 16 months ago by
Replying to egourgoulhon:
Thanks for improving significantly the generator of nonredundant 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 followup: 11 Changed 16 months ago by
Replying to ghspaghettisalat:
Replying to egourgoulhon:
Thanks for improving significantly the generator of nonredundant 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 Changed 16 months ago by
Replying to egourgoulhon:
Replying to ghspaghettisalat:
Replying to egourgoulhon:
Thanks for improving significantly the generator of nonredundant 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
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.
comment:13 Changed 16 months ago by
Commit:  1c108d85c05ccd171548f14f21b1c5a7db896d93 → 08ccd3b89d25893230bbef55b0ed448b23c9d3de 

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
08ccd3b  tensor module components: nonredundant index generator: performance improvements

comment:14 Changed 16 months ago by
I have implemented sorting of the indices for display. All outputs should now appear the same way as before.
Changed 6 hours ago by ghmjungmath:
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 nonsymmetrized 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
Status:  needs_work → needs_review 

comment:16 Changed 16 months ago by
Reviewers:  → Eric Gourgoulhon, Michael Jung 

Status:  needs_review → positive_review 
Thanks for the changes! Michael, do you agree with the positive review?
comment:17 followup: 19 Changed 16 months ago by
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:19 followup: 22 Changed 16 months ago by
Replying to ghmjungmath:
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 followup: 23 Changed 16 months ago by
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 nonsymmetrized indices as being symmetrized with themselves + # treat nonsymmetrized 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
I ponder, is the outer whileloop 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.
comment:22 Changed 16 months ago by
Replying to ghspaghettisalat:
Replying to ghmjungmath:
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:23 Changed 16 months ago by
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
Milestone:  sage9.4 → sage9.5 

comment:25 Changed 15 months ago by
Branch:  public/tensorindices → 08ccd3b89d25893230bbef55b0ed448b23c9d3de 

Resolution:  → fixed 
Status:  positive_review → closed 
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
tensor module components: nonredundant index generator: performance improvements