Opened 3 months ago

Closed 7 weeks ago

#32318 closed enhancement (fixed)

Generate non-redundant indices of tensors with symmetries efficiently

Reported by: gh-spaghettisalat Owned by:
Priority: major Milestone: sage-9.5
Component: performance Keywords:
Cc: slelievre, gh-mjungmath, 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:

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 3 months ago by gh-spaghettisalat

  • Status changed from new to needs_review

comment:2 Changed 3 months ago by gh-spaghettisalat

  • Summary changed from Performace improvements in generating non-redundant indices of tensors with indices to Performace improvements in generating non-redundant indices of tensors with symmetries

comment:3 Changed 3 months ago by slelievre

  • Cc slelievre added
  • Summary changed from Performace improvements in generating non-redundant indices of tensors with symmetries to Generate non-redundant indices of tensors with symmetries efficiently

comment:4 Changed 3 months ago by git

  • Commit changed from a60d93c37102245cb7f9e319576d9f92e7d15c1e to 1c108d85c05ccd171548f14f21b1c5a7db896d93

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 3 months ago by mkoeppe

  • Cc gh-mjungmath egourgoulhon added

comment:6 Changed 3 months 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.
  • Very minor (PEP norm):
    - for k in range(1,len(isym)):
    + for k in range(1, len(isym)):
    

comment:7 Changed 3 months ago by gh-mjungmath

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 3 months ago by gh-mjungmath (previous) (diff)

comment:8 follow-up: Changed 3 months ago by egourgoulhon

  • Status changed from needs_review to needs_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 ; follow-up: Changed 3 months ago by 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? 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 ; follow-up: Changed 3 months ago by 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 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 3 months ago by gh-spaghettisalat

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 3 months ago by gh-mjungmath

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

Last edited 3 months ago by gh-mjungmath (previous) (diff)

comment:13 Changed 3 months ago by git

  • Commit changed from 1c108d85c05ccd171548f14f21b1c5a7db896d93 to 08ccd3b89d25893230bbef55b0ed448b23c9d3de

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 3 months ago by gh-spaghettisalat

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 3 months ago by gh-spaghettisalat

  • Status changed from needs_work to needs_review

comment:16 Changed 3 months ago by egourgoulhon

  • 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 follow-up: Changed 3 months ago by 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.

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

comment:18 Changed 3 months ago by gh-mjungmath

This is a huge improvement though. Thank you Marius!

comment:19 in reply to: ↑ 17 ; follow-up: Changed 3 months ago by 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...

comment:20 follow-up: Changed 3 months ago by slelievre

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 3 months ago by gh-mjungmath

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 3 months ago by gh-mjungmath (previous) (diff)

comment:22 in reply to: ↑ 19 Changed 3 months ago by gh-mjungmath

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 3 months ago by egourgoulhon

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 3 months ago by mkoeppe

  • Milestone changed from sage-9.4 to sage-9.5

comment:25 Changed 7 weeks ago by vbraun

  • Branch changed from public/tensor-indices to 08ccd3b89d25893230bbef55b0ed448b23c9d3de
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.