Tensor contraction in SymPy 1.9 take too long

[K4stan]

After i upgrade from SymPy 1.4 -> 1.9 The script that took about 30 min to finish on my PC now take like forever. The problem is tensorcontraction function, it seems it take much longer than with SymPy 1.4. IS there some known change that can cause this problem ?

[moorepants]

Please provide a minimal example that demonstrates the issue. That would allow sympy developers to investigate.

[oscarbenjamin]

From the looks of things this is converting a tuple into an integer and then parsing it back into a tuple:

sympy/sympy/tensor/array/arrayop.py

Lines 181 to 186 in 666c4f9

	for icontrib in itertools.product(*remaining_indices):
	index_base_position = sum(icontrib)
	isum = S.Zero
	for sum_to_index in itertools.product(*summed_deltas):
	idx = array._get_tuple_index(index_base_position + sum(sum_to_index))
	isum += array[idx]

[oscarbenjamin]

And then __getitem__ calls _parse_index to convert the tuple back into an int so that it can be used to index into the internal flat list. Bypassing all of these speeds things up massively:

diff --git a/sympy/tensor/array/arrayop.py b/sympy/tensor/array/arrayop.py
index c82b0d8..a83b7fc 100644
--- a/sympy/tensor/array/arrayop.py
+++ b/sympy/tensor/array/arrayop.py
@@ -182,8 +182,8 @@ def tensorcontraction(array, *contraction_axes):
         index_base_position = sum(icontrib)
         isum = S.Zero
         for sum_to_index in itertools.product(*summed_deltas):
-            idx = array._get_tuple_index(index_base_position + sum(sum_to_index))
-            isum += array[idx]
+            idx = index_base_position + sum(sum_to_index)
+            isum += array._array[idx]
 
         contracted_array.append(isum)

[Upabjojr]

@K4stan I was looking at your example:

T_Gr=syp.tensorproduct(C2x,C2x,C2x,C2x,T_G)
for ii in range(4):
    T_Gr=syp.tensorcontraction(T_Gr,(ii+1,8-ii))

I can suggest to use as a temporary solution to use a lazy version of the array contraction:

In [2]: from sympy.tensor.array.expressions.array_expressions import *

In [3]: T_Gr = ArrayTensorProduct(C2x,C2x,C2x,C2x,T_G)

In [4]: for ii in range(4):
    ...:     T_Gr = ArrayContraction(T_Gr, (ii+1, 8-ii)).doit()
    ...:

In [5]: T_Gr
Out[5]: ArrayContraction(ArrayTensorProduct([[[[G1111, G1112, G1113], [G1112, G1122, G1123], [G1113, G1123, G1133]], [[G1211, G1212, G1213], [G1212, G1222, G1223], [G1213, G1223, G1233]], [[G1311, G1312, G1313], [G1312, G1322, G1323], [G1313, G1323, G1333]]], [[[G1211, G1212, G1213], [G1212, G1222, G1223], [G1213, G1223, G1233]], [[G2211, G2212, G2213], [G2212, G2222, G2223], [G2213, G2223, G2233]], [[G2311, G2312, G2313], [G2312, G2322, G2323], [G2313, G2323, G2333]]], [[[G1311, G1312, G1313], [G1312, G1322, G1323], [G1313, G1323, G1333]], [[G2311, G2312, G2313], [G2312, G2322, G2323], [G2313, G2323, G2333]], [[G3311, G3312, G3313], [G3312, G3322, G3323], [G3313, G3323, G3333]]]], [[1, 0, 0], [0, -1, 0], [0, 0, -1]], [[1, 0, 0], [0, -1, 0], [0, 0, -1]], [[1, 0, 0], [0, -1, 0], [0, 0, -1]], [[1, 0, 0], [0, -1, 0], [0, 0, -1]]), (0, 5), (1, 7), (2, 9), (3, 11))

In [6]: T_Gr.as_explicit()
Out[6]: [[[[G1111, -G1112, -G1113], [-G1112, G1122, G1123], [-G1113, G1123, G1133]], [[-G1211, G1212, G1213], [G1212, -G1222, -G1223], [G1213, -G1223, -G1233]], [[-G1311, G1312, G1313], [G1312, -G1322, -G1323], [G1313, -G1323, -G1333]]], [[[-G1211, G1212, G1213], [G1212, -G1222, -G1223], [G1213, -G1223, -G1233]], [[G2211, -G2212, -G2213], [-G2212, G2222, G2223], [-G2213, G2223, G2233]], [[G2311, -G2312, -G2313], [-G2312, G2322, G2323], [-G2313, G2323, G2333]]], [[[-G1311, G1312, G1313], [G1312, -G1322, -G1323], [G1313, -G1323, -G1333]], [[G2311, -G2312, -G2313], [-G2312, G2322, G2323], [-G2313, G2323, G2333]], [[G3311, -G3312, -G3313], [-G3312, G3322, G3323], [-G3313, G3323, G3333]]]]

In this code example, ArrayContraction creates a lazy version of tensorcontraction (that is, it creates the contraction structure, but does not evaluate). Calling .doit() on ArrayContraction objects will flatten nested ArrayContraction objects into a single one, correcting for the contraction indices.

At this point, calling .as_explicit() on the final object will call the tensor contraction algorithm only once.

You can see the new contraction indices by calling:

In [7]: T_Gr.contraction_indices
Out[7]: ((0, 5), (1, 7), (2, 9), (3, 11))

[K4stan]

Ok, I see, thanks for your input.

[Upabjojr]

Even better would be to integrate the opt_einsum library into SymPy.