# Perform matrix/tensor contractions more efficiently

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?

• Are you absolutely certain that they can't be done using Dot (which is more general than just matrix multiplication), Transpose, Flatten, and Tr? I'm willing to bet that it can be done with those optimized built-ins. Also, at least, please provide an example of one that can't be done with matrix multiplication and Tr for us to play with. – march Mar 30 '16 at 21:12
• The important fact you have to know is that Dot[] always contracts the last index of one tensor into the first of another. Therefore, you just need to Transpose or Flatten appropriately and then use Dot. As @march said this will always work and you can do the most complicated contractions. – Michael Weyrauch Mar 30 '16 at 22:27
• You've seen MatrixPower[], right? – J. M.'s torpor Mar 30 '16 at 22:58

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

• Thanks for the answer, I actually only now really got to use it! Could you maybe explain why I get an out of memory message if I don't use Inactive for the TensorProduct? – Stan May 19 '16 at 23:37
• The tensor product of small arrays can be a very large array. The tensor contraction will then reduce that large array to a small result, and the use of Inactive allows getting to that final result without actually constructing the intermediate large array. If you don't use Inactive, TensorProduct will construct a huge array, and that's why you get the memory error. In my example above A has 100^3 = 10^6 entries and A4 would have 100^12 = 10^24 entries if allowed to expand, which is too much for computer memory. – jose May 24 '16 at 3:18

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 :)