Opened 10 months ago
Closed 8 months ago
#32318 closed enhancement (fixed)
Generate nonredundant indices of tensors with symmetries efficiently
Reported by:  ghspaghettisalat  Owned by:  

Priority:  major  Milestone:  sage9.5 
Component:  performance  Keywords:  
Cc:  slelievre, ghmjungmath, egourgoulhon, tscrim  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 10 months ago by
 Status changed from new to needs_review
comment:2 Changed 10 months ago by
 Summary changed from Performace improvements in generating nonredundant indices of tensors with indices to Performace improvements in generating nonredundant indices of tensors with symmetries
comment:3 Changed 10 months ago by
 Cc slelievre added
 Summary changed from Performace improvements in generating nonredundant indices of tensors with symmetries to Generate nonredundant indices of tensors with symmetries efficiently
comment:4 Changed 10 months ago by
 Commit changed from a60d93c37102245cb7f9e319576d9f92e7d15c1e to 1c108d85c05ccd171548f14f21b1c5a7db896d93
comment:5 Changed 10 months ago by
 Cc ghmjungmath egourgoulhon added
comment:6 Changed 10 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 10 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 10 months ago by
 Status changed from needs_review to 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 in reply to: ↑ 8 ; followup: ↓ 10 Changed 10 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 in reply to: ↑ 9 ; followup: ↓ 11 Changed 10 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 in reply to: ↑ 10 Changed 10 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 10 months ago by
 Cc tscrim 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 10 months ago by
 Commit changed from 1c108d85c05ccd171548f14f21b1c5a7db896d93 to 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 10 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 10 months ago by
 Status changed from needs_work to needs_review
comment:16 Changed 9 months ago by
 Reviewers set to Eric Gourgoulhon, Michael Jung
 Status changed from needs_review to positive_review
Thanks for the changes! Michael, do you agree with the positive review?
comment:17 followup: ↓ 19 Changed 9 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:18 Changed 9 months ago by
This is a huge improvement though. Thank you Marius!
comment:19 in reply to: ↑ 17 ; followup: ↓ 22 Changed 9 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 9 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 9 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 in reply to: ↑ 19 Changed 9 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 in reply to: ↑ 20 Changed 9 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 9 months ago by
 Milestone changed from sage9.4 to sage9.5
comment:25 Changed 8 months ago by
 Branch changed from public/tensorindices to 08ccd3b89d25893230bbef55b0ed448b23c9d3de
 Resolution set to fixed
 Status changed from positive_review to 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