Bug in multiple tensor contractions with scalar result

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

CC: @tscrim @mjungmath @mkoeppe

Component: manifolds

Keywords: tensor contraction

Author: Eric Gourgoulhon

Branch/Commit: 48f5523

Reviewer: Michael Jung

Issue created by migration from https://trac.sagemath.org/ticket/32355

[egourgoulhon]

comment:2

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[[...]].

[egourgoulhon]

comment:3

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)

[egourgoulhon]

Branch: public/manifolds/contraction_bug-32355

[egourgoulhon]

Commit: 26ba8f7

[egourgoulhon]

Author: Eric Gourgoulhon

[mjungmath]

comment:5

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.

[sagetrac-git]

Changed commit from 26ba8f7 to 3e5ac7d

[sagetrac-git]

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

3e5ac7d

Improvement in the fix for #32355

[egourgoulhon]

comment:7

Replying to @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...

[mjungmath]

comment:8

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.

[sagetrac-git]

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

48f5523

Yet another improvement in the fix for #32355

[sagetrac-git]

Changed commit from 3e5ac7d to 48f5523

[egourgoulhon]

comment:10

Replying to @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.

[mjungmath]

Reviewer: Michael Jung

[egourgoulhon]

comment:12

Thank you for the review!

[vbraun]

Changed branch from public/manifolds/contraction_bug-32355 to 48f5523