Compiling einstein sums optimally Einstein summation is a convenient way to express tensor operations which has found its way in tensor libraries like numpy, torch, tensorflow, etc.
Its flexibility lets us represent the product of three matrices, $X$, $Y$, $Z$  of dimensions $(a,b)$, $(b,c)$, $(c,d)$ as
X.Y.Z = einsum('ab,bc,cd-&gt;ad',X,Y,Z)

However, the above compiles to something like
for a_ in range(a):
  for d_ in range(d):
     res[a_,d_] = 0
     for b_ in range(b):
        for c_ in range(c):
           res[a_,d_] += X[a_,b_] * Y[b_,c_] * Z[c_, d_]

This native is quadratic in the size of the matrices when simply doing
einsum('ac,cd-&gt;ad',einsum('ab,bc'-&gt;'ac', X, Y), Z)

Would be merely cubic.
There are roughly three levels of optimization we can imagine a smarter implementation of einsum to perform.

Decompose an einsum of tensors $(x_1, \ldots, x_n)$ into an einsum of pairs of tensors $x_1, x_2$, $e(x_1, x_2), x_3$, etc to optimize computation time.

Rely on associativity (where applicable) to pick those pairs judiciously (this is a classic dynamic programming problem) and construct the appropriate intermediary tensors.

Discover Strassen-like formulas for the specific tensor computation


While 3 seems clearly out of reach, 1 and 2 seem like they could be achievable exactly with a reasonably straightforward algorithm. Are such algorithms known for generic einstein summations? Have they been studied?
It seems that the general problem of finding the optimal contraction order is NP-hard [1]. A recent paper on approximately optimizing the contraction order, and containing relevant references, is [2].
[1] Chi-Chung, Lam, P. Sadayappan, and Rephael Wenger. &quot;On optimizing a class of multi-dimensional loops with reduction for parallel execution.&quot; Parallel Processing Letters 7.02 (1997): 157-168.
[2] Schindler, Frank, and Adam Jermyn. &quot;Algorithms for tensor network contraction ordering.&quot; Machine Learning: Science and Technology (2020).