<p>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].</p>
<p>[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.</p>
<p>[2] Schindler, Frank, and Adam Jermyn. &quot;Algorithms for tensor network contraction ordering.&quot; Machine Learning: Science and Technology (2020).</p>

CC BY-SA 4.0
2020-08-20T10:30:14.167Z
Thank you, "tensor network contraction" is the keyword I was missing
Arthur B
CC BY-SA 4.0
2020-08-20T10:38:11.077Z
Not all evaluation strategies are contractions though. The naive product of 3 $n \times n$ matrices in $\mathcal{O}(n^4)$ time doesn't happen if you follow *any* one of the two contraction strategy... yet, you'll see numpy's einsum for instance run it in quadratic time anyway. What gives?
Arthur B
CC BY-SA 4.0
2020-08-20T10:49:55.973Z
By quadratic, do you actually mean *quartic* (${\cal O}(n^4)$)? I guess that any *optimal* evaluation strategies would consist of contractions.
smapers
CC BY-SA 4.0
2020-08-20T17:39:45.127Z
Sorry, yes I meant to write quartic. And the answer to my question is that, yes, this is a thing, but there's a flag...
https://optimized-einsum.readthedocs.io/en/stable/index.html
Arthur B
CC BY-SA 4.0
2020-08-20T17:41:33.370Z
Specifically I was thinking: the paper assumes that contractions always dominate naive evaluation, but surely if that were the case, numpy's einsum would implement at least a naive contraction. Turn out it does, but only with a flag.
Arthur B
CC BY-SA 4.0
2020-08-19T09:33:50.377Z
2020-08-19T09:33:50.377Z
smapers
<p>Einstein summation is a convenient way to express tensor operations which has found its way in tensor libraries like numpy, torch, tensorflow, etc.</p>
<p>Its flexibility lets us represent the product of three matrices, <span class="math-container">$X$</span>, <span class="math-container">$Y$</span>, <span class="math-container">$Z$</span>  of dimensions <span class="math-container">$(a,b)$</span>, <span class="math-container">$(b,c)$</span>, <span class="math-container">$(c,d)$</span> as</p>
<pre><code>X.Y.Z = einsum('ab,bc,cd-&gt;ad',X,Y,Z)
</code></pre>
<p>However, the above compiles to something like</p>
<pre><code>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_]
</code></pre>
<p>This native is quadratic in the size of the matrices when simply doing</p>
<pre><code>einsum('ac,cd-&gt;ad',einsum('ab,bc'-&gt;'ac', X, Y), Z)
</code></pre>
<p>Would be merely cubic.</p>
<p>There are roughly three levels of optimization we can imagine a smarter implementation of <code>einsum</code> to perform.</p>
<ol>
<li><p>Decompose an einsum of tensors <span class="math-container">$(x_1, \ldots, x_n)$</span> into an einsum of pairs of tensors <span class="math-container">$x_1, x_2$</span>, <span class="math-container">$e(x_1, x_2), x_3$</span>, etc to optimize computation time.</p>
</li>
<li><p>Rely on associativity (where applicable) to pick those pairs judiciously (this is a classic dynamic programming problem) and construct the appropriate intermediary tensors.</p>
</li>
<li><p>Discover Strassen-like formulas for the specific tensor computation</p>
</li>
</ol>
<p>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?</p>

CC BY-SA 4.0
2020-08-18T15:53:12.847Z
What does the notation `cd->ab` mean?
Mark Schultz-Wu
CC BY-SA 4.0
2020-08-18T21:36:40.197Z
The quote inserted was a typo, it's `ab,cd->ad`
Arthur B
<optimization><matrices>
Compiling einstein sums optimally
Arthur B
# Compiling einstein sums optimally

(Asked by: username_0 on 2020-08-18)

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->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->ad',einsum('ab,bc'->'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?

**Tags:** optimization, matrices

> **Comments:**
> **username_1**: What does the notation `cd->ab` mean?  
> **username_0**: The quote inserted was a typo, it's `ab,cd->ad`  

---

## 1 answers

### Answer by username_2 on 2020-08-19. Score: 2 | Accepted answer

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. "On optimizing a class of multi-dimensional loops with reduction for parallel execution." Parallel Processing Letters 7.02 (1997): 157-168.

[2] Schindler, Frank, and Adam Jermyn. "Algorithms for tensor network contraction ordering." Machine Learning: Science and Technology (2020).

> **Comments:**
> **username_0**: Thank you, "tensor network contraction" is the keyword I was missing  
> **username_0**: Not all evaluation strategies are contractions though. The naive product of 3 $n \times n$ matrices in $\mathcal{O}(n^4)$ time doesn't happen if you follow *any* one of the two contraction strategy... yet, you'll see numpy's einsum for instance run it in quadratic time anyway. What gives?  
> **username_2**: By quadratic, do you actually mean *quartic* (${\cal O}(n^4)$)? I guess that any *optimal* evaluation strategies would consist of contractions.  
> **username_0**: Sorry, yes I meant to write quartic. And the answer to my question is that, yes, this is a thing, but there's a flag... https://optimized-einsum.readthedocs.io/en/stable/index.html  
> **username_0**: Specifically I was thinking: the paper assumes that contractions always dominate naive evaluation, but surely if that were the case, numpy's einsum would implement at least a naive contraction. Turn out it does, but only with a flag.  

---

