Opened 3 months ago

Closed 7 weeks ago

#32355 closed defect (fixed)

Bug in multiple tensor contractions with scalar result

Reported by: egourgoulhon Owned by:
Priority: major Milestone: sage-9.5
Component: manifolds Keywords: tensor contraction
Cc: tscrim, gh-mjungmath, mkoeppe Merged in:
Authors: Eric Gourgoulhon Reviewers: Michael Jung
Report Upstream: N/A Work issues:
Branch: 48f5523 (Commits, GitHub, GitLab) Commit: 48f5523074fbbb835fa41be1d567d244e9244ae4
Dependencies: Stopgaps:

Status badges

Description (last modified by egourgoulhon)

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

  • Description modified (diff)
  • Summary changed from Bug in multiple tensor contractions to Bug in multiple tensor contractions with scalar result

comment:2 Changed 3 months ago by egourgoulhon

The culprit is line 2187 of src/sage/tensor/modules/comp.py:

    for ind in comp_for_contr.index_generator():
        res += self[[ind]] * other[[ind]]

ind should be reordered before being passed to other[[...]].

comment:3 Changed 2 months ago by egourgoulhon

  • Branch set to public/manifolds/contraction_bug-32355
  • Cc tscrim gh-mjungmath 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:

26ba8f7Fix bug in Components.contract() (Trac #32355)

comment:4 Changed 2 months ago by egourgoulhon

  • Authors set to Eric Gourgoulhon

comment:5 follow-up: Changed 2 months ago by gh-mjungmath

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 for-loop into a list comprehension which would be faster but I reckon less well readable.

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

comment:6 Changed 2 months ago by git

  • Commit changed from 26ba8f76d35623c499ce7c70a0b6a2afbb724929 to 3e5ac7d48bc15509694e0b904b8595b38caa4e8a

Branch pushed to git repo; I updated commit sha1. New commits:

3e5ac7dImprovement in the fix for #32355

comment:7 in reply to: ↑ 5 ; follow-up: Changed 2 months ago by egourgoulhon

Replying to gh-mjungmath:

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 ; follow-up: Changed 2 months ago by gh-mjungmath

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.

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

comment:9 Changed 2 months ago by git

  • Commit changed from 3e5ac7d48bc15509694e0b904b8595b38caa4e8a to 48f5523074fbbb835fa41be1d567d244e9244ae4

Branch pushed to git repo; I updated commit sha1. New commits:

48f5523Yet another improvement in the fix for #32355

comment:10 in reply to: ↑ 8 Changed 2 months ago by egourgoulhon

Replying to gh-mjungmath: 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 2 months ago by gh-mjungmath

  • Reviewers set to Michael Jung
  • Status changed from needs_review to positive_review

comment:12 Changed 2 months ago by egourgoulhon

Thank you for the review!

comment:13 Changed 2 months ago by mkoeppe

  • Milestone changed from sage-9.4 to sage-9.5

comment:14 Changed 7 weeks ago by vbraun

  • Branch changed from public/manifolds/contraction_bug-32355 to 48f5523074fbbb835fa41be1d567d244e9244ae4
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.