Bug in multiple tensor contractions with scalar result

Reported by: Owned by: Eric Gourgoulhon major sage-9.5 manifolds tensor contraction Travis Scrimshaw, Michael Jung, Matthias Köppe Eric Gourgoulhon Michael Jung N/A 48f5523 48f5523074fbbb835fa41be1d567d244e9244ae4

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.

comment:1 Changed 16 months ago by Eric Gourgoulhon

Description: modified (diff) Bug in multiple tensor contractions → Bug in multiple tensor contractions with scalar result

comment:2 Changed 16 months ago by Eric Gourgoulhon

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 16 months ago by Eric Gourgoulhon

Branch: → public/manifolds/contraction_bug-32355 Travis Scrimshaw Michael Jung Matthias Köppe added → 26ba8f76d35623c499ce7c70a0b6a2afbb724929 new → 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 16 months ago by Eric Gourgoulhon

Authors: → Eric Gourgoulhon

comment:5 follow-up:  7 Changed 16 months ago by Michael Jung

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

```-            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 16 months ago by Michael Jung (previous) (diff)

comment:6 Changed 16 months ago by git

Commit: 26ba8f76d35623c499ce7c70a0b6a2afbb724929 → 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 ; follow-up:  8 Changed 16 months ago by Eric Gourgoulhon

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.

```-            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:  10 Changed 16 months ago by Michael Jung

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 16 months ago by Michael Jung (previous) (diff)

comment:9 Changed 16 months ago by git

Commit: 3e5ac7d48bc15509694e0b904b8595b38caa4e8a → 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 16 months ago by Eric Gourgoulhon

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 16 months ago by Michael Jung

Reviewers: → Michael Jung needs_review → positive_review

comment:12 Changed 16 months ago by Eric Gourgoulhon

Thank you for the review!

comment:13 Changed 16 months ago by Matthias Köppe

Milestone: sage-9.4 → sage-9.5

comment:14 Changed 15 months ago by Volker Braun

Branch: public/manifolds/contraction_bug-32355 → 48f5523074fbbb835fa41be1d567d244e9244ae4 → fixed positive_review → closed
Note: See TracTickets for help on using tickets.