Opened 9 months ago
Closed 8 months ago
#32355 closed defect (fixed)
Bug in multiple tensor contractions with scalar result
Reported by:  egourgoulhon  Owned by:  

Priority:  major  Milestone:  sage9.5 
Component:  manifolds  Keywords:  tensor contraction 
Cc:  tscrim, ghmjungmath, mkoeppe  Merged in:  
Authors:  Eric Gourgoulhon  Reviewers:  Michael Jung 
Report Upstream:  N/A  Work issues:  
Branch:  48f5523 (Commits, GitHub, GitLab)  Commit:  48f5523074fbbb835fa41be1d567d244e9244ae4 
Dependencies:  Stopgaps: 
Description (last modified by )
Tensor contractions on multiple indices give wrong results when the outcome is a scalar and the index positions in each tensor do not have the same order. For instance, let us consider
sage: M = Manifold(2, 'M') sage: X.<x,y> = M.chart() sage: a = M.multivector_field(2) sage: a[0,1] = 1 sage: b = M.diff_form(2) sage: b[0,1] = 1
The double contraction a^{ij} b_{ij} is obtained via
sage: a.contract(0, 1, b, 0, 1).expr() 2 sage: (a['^ij']*b['_ij']).expr() # equivalent to above 2
which is correct. However, asking for the double contraction a^{ij} b_{ji} leads to
sage: a.contract(0, 1, b, 1, 0).expr() 2 sage: (a['^ij']*b['_ji']).expr() # equivalent to above 2
which is obviously wrong: the result should be 2, since b is antisymmetric.
If the outcome is not a scalar field, there is no issue:
sage: c = M.tensor_field(3, 0, antisym=(0,1)) sage: c[0,1,0] = 1 sage: c.contract(0, 1, b, 0, 1).display() # correct result 2 ∂/∂x sage: c.contract(0, 1, b, 1, 0).display() # correct result 2 ∂/∂x
This bug has been reported in https://ask.sagemath.org/question/58379.
Change History (14)
comment:1 Changed 9 months ago by
 Description modified (diff)
 Summary changed from Bug in multiple tensor contractions to Bug in multiple tensor contractions with scalar result
comment:2 Changed 9 months ago by
comment:3 Changed 9 months ago by
 Branch set to public/manifolds/contraction_bug32355
 Cc tscrim ghmjungmath mkoeppe added
 Commit set to 26ba8f76d35623c499ce7c70a0b6a2afbb724929
 Status changed from new to needs_review
The bug is fixed in the above branch. Some slight clean up of the code for parallelized contractions has also been performed.
New commits:
26ba8f7  Fix bug in Components.contract() (Trac #32355)

comment:4 Changed 9 months ago by
comment:5 followup: ↓ 7 Changed 9 months ago by
Here are some comments.
 I think you don't need to initialize a new component in
comp_for_contr = Components(self._ring, self._frame, ncontr, start_index=self._sindex)
and you can use
self.index_generator()
right away, no?
 What about this?
 for ind_s in comp_for_contr.index_generator():  ind_o = [None for i in range(ncontr)]  for pos_s, pos_o in contractions:  ind_o[pos_o] = ind_s[pos_s]  ind_pairs.append((ind_s, ind_o)) + for ind_s in self.index_generator(): + ind_o = [ind_s[pos_s] for pos_s, _ in sorted(contractions, key=itemgetter(1))] + ind_pairs.append((ind_s, ind_o))
You could also turn the outer forloop into a list comprehension which would be faster but I reckon less well readable.
comment:6 Changed 9 months ago by
 Commit changed from 26ba8f76d35623c499ce7c70a0b6a2afbb724929 to 3e5ac7d48bc15509694e0b904b8595b38caa4e8a
Branch pushed to git repo; I updated commit sha1. New commits:
3e5ac7d  Improvement in the fix for #32355

comment:7 in reply to: ↑ 5 ; followup: ↓ 8 Changed 9 months ago by
Replying to ghmjungmath:
Here are some comments.
Thanks for taking a look into this.
 I think you don't need to initialize a new component in
comp_for_contr = Components(self._ring, self._frame, ncontr, start_index=self._sindex)and you can use
self.index_generator()
right away, no?
Yes indeed! (this is needed for the generic case, but not for the specific case of a scalar output). The latest commit implement this change.
 What about this?
 for ind_s in comp_for_contr.index_generator():  ind_o = [None for i in range(ncontr)]  for pos_s, pos_o in contractions:  ind_o[pos_o] = ind_s[pos_s]  ind_pairs.append((ind_s, ind_o)) + for ind_s in self.index_generator(): + ind_o = [ind_s[pos_s] for pos_s, _ in sorted(contractions, key=itemgetter(1))] + ind_pairs.append((ind_s, ind_o))
I would prefer to keep the original version, since it is more efficient from a computational point of view: there is no need to sort a list. Granted: in most use cases, the list contains only a few elements and sorting shall be fast. So it is more a matter of taste...
comment:8 in reply to: ↑ 7 ; followup: ↓ 10 Changed 9 months ago by
Replying to egourgoulhon:
I would prefer to keep the original version, since it is more efficient from a computational point of view: there is no need to sort a list. Granted: in most use cases, the list contains only a few elements and sorting shall be fast. So it is more a matter of taste...
You're right, this is faster! Then I still have one last comment:
 ind_o = [None for i in range(ncontr)] + ind_o = [None] * ncontr
This is indeed faster.
If that is done, and patchbot is happy > positive review.
comment:9 Changed 9 months ago by
 Commit changed from 3e5ac7d48bc15509694e0b904b8595b38caa4e8a to 48f5523074fbbb835fa41be1d567d244e9244ae4
Branch pushed to git repo; I updated commit sha1. New commits:
48f5523  Yet another improvement in the fix for #32355

comment:10 in reply to: ↑ 8 Changed 9 months ago by
Replying to ghmjungmath: Then I still have one last comment:
 ind_o = [None for i in range(ncontr)] + ind_o = [None] * ncontr
Done in the above commit.
comment:11 Changed 9 months ago by
 Reviewers set to Michael Jung
 Status changed from needs_review to positive_review
comment:12 Changed 9 months ago by
Thank you for the review!
comment:13 Changed 9 months ago by
 Milestone changed from sage9.4 to sage9.5
comment:14 Changed 8 months ago by
 Branch changed from public/manifolds/contraction_bug32355 to 48f5523074fbbb835fa41be1d567d244e9244ae4
 Resolution set to fixed
 Status changed from positive_review to closed
The culprit is line 2187 of
src/sage/tensor/modules/comp.py
:ind
should be reordered before being passed toother[[...]]
.