Perform matrix/tensor contractions more efficiently

Question

Problem: explicit index contractions for matrices using the Sum command take too long, and I want to improve the performance for complicated computations.

Let me give you an example, let M be a random $8\times 8$ matrix:

M = RandomReal[{0, 1}, {8, 8}];

If I want to multiply 8 of those matrices and take a trace I would write something like

Tr[M.M.M.M.M.M.M.M]

This gives me the result almost instantaneously.

I can write the same computation in a more complicated way by using the Sum or equivalently ParallelSum command and writing out the index contractions explicitly:

ParallelSum[
M[[i1, i2]] M[[i2, i3]] M[[i3, i4]] M[[i4, i5]]
M[[i5, i6]] M[[i6, i7]] M[[i7, i8]] M[[i8, i1]],
{i1, 1, 8}, {i2, 1, 8}, {i3, 1, 8}, {i4, 1, 8},
{i5, 1, 8}, {i6, 1, 8}, {i7, 1, 8}, {i8, 1, 8}]

This computation takes literally ages, and I actually never waited till the end to see the final result.

Question: how can I implement explicit contractions like in the last example efficiently? I need to do this because in more complicated computations with multidimensional tensors contractions cannot be written using matrix multiplication and traces. What is the best practise for this kind of computations?

Answer 1

I think the simplest way to handle the general case is to use TensorProduct and TensorContract, as follows:

Take a rank 3 array for example, in dimension 100:

In[1]:= A = RandomReal[{-1, 1}, {100, 100, 100}];

Construct a rank 12 array. Note the use of Inactive, to avoid TensorProduct constructing a large intermediate array:

In[2]:= A4 = Inactive[TensorProduct][A, A, A, A];

In[3]:= TensorRank[A4]
Out[3]= 12

Now choose any random contraction of all levels:

In[4]:= contraction = Partition[RandomSample[Range[12]], 2]
Out[4]= {{1, 3}, {11, 6}, {8, 7}, {12, 9}, {4, 5}, {2, 10}}

That means a contraction of the 1st and 3rd levels, the 11th and 6th levels, etc. You can easily see the levels as indices.

Finally TensorContract will perform the contraction, and quite fast. Note again the use of Activate to eliminate the Inactive head:

In[5]:= TensorContract[A4, contraction] // Activate // AbsoluteTiming
Out[5]= {0.006712, -115212.}

Answer 2

I am illustrating my answer with a very short code. Rather, it is a long comment. The computational complexity of your two examples is cardinally different. In the first case you perform 7 matrix multiplications and perform a trace on the result. The last operation has a quadratic scaling. Therefore leading complexity is $7 N^3$, where $N$ is the matrix dimension.

In your second example you do all multiplications at once. Computational complexity is $N^8$. Quite clear $N^8\gg 7N^3$ !

Now, can we do even better? Yes it is indeed possible. You can do the following sequence of operations:

 M=M.M
 M=M.M
 M=M.M
 Tr[M]

The computational complexity is $3 N^3$. We have a speed up of 7/3.

Can we do even better? Yes it is indeed possible. In the above example $3\times 8^3=1536$ operations (count only multiplications) is required. We can use the simplest fast matrix multiplication algorithm due to Strassen. It requires 7 multiplications to multiply matrix $2\times2$. Matrices with dimensions that are powers of 2 are well suitable to apply the algorithm recursively. To multiply 2 matrices $8\times8$ only $7^3=343$ operations are required. In total we need only 1029 multiplications. Considering that we only need diagonal elements for the trace we can reduce the estimate even further $2\times 7^3+8^2=750$.

Can we do even better? Quite probably yes, however, the algorithm for fast matrix multiplications of matrices $4\times4$ with less then 49 multiplications is currently unknown. It was conjectured, however, that asymptotically ($N\rightarrow \infty$) the scaling is quadratic.

To summarise:

$16777216\gg1536>1029>750\ge$computational complexity$>192$. Thus, there is still a lot of room for improvement :)