# Large Expressions with Greedy

Using the greedy method allows the contraction of hundreds of tensors. Here's an example from quantum of computing the inner product between two 'Matrix Product States'. Graphically, if we represent each tensor as an O, give it the same number of 'legs' as it has indices, and join those legs when that index is summed with another tensor, we get an expression for n particles that looks like:

O-O-O-O-O-O-     -O-O-O-O-O-O
| | | | | |  ...  | | | | | |
O-O-O-O-O-O-     -O-O-O-O-O-O

0 1 2 3 4 5 ........... n-2 n-1


The meaning of this is not that important other than its a large, useful contraction. For n=100 it involves 200 different tensors and about 300 unique indices. With this many indices it can be useful to generate them with the function opt_einsum.parser.get_symbol.

### Setup the string

import numpy as np
import opt_einsum as oe

n = 100
phys_dim = 3
bond_dim = 10

# O--
# |
# O--
einsum_str = "ab,ac,"

for i in range(1, n - 1):
# set the upper left/right, middle and lower left/right indices
# --O--
#   |
# --O--
j = 3 * i
ul, ur, m, ll, lr = (oe.get_symbol(i)
for i in (j - 1, j + 2, j, j - 2, j + 1))
einsum_str += "{}{}{},{}{}{},".format(m, ul, ur, m, ll, lr)

# finish with last site
# --O
#   |
# --O
i = n - 1
j = 3 * i
ul, m, ll, =  (oe.get_symbol(i) for i in (j - 1, j, j - 2))
einsum_str += "{}{},{}{}".format(m, ul, m, ll)


### Generate the shapes

def gen_shapes():
yield (phys_dim, bond_dim)
yield (phys_dim, bond_dim)
for i in range(1, n - 1):
yield(phys_dim, bond_dim, bond_dim)
yield(phys_dim, bond_dim, bond_dim)
yield (phys_dim, bond_dim)
yield (phys_dim, bond_dim)

shapes = tuple(gen_shapes())


Let's time how long it takes to generate the expression ('greedy' is used by default, and we turn off the memory_limit):

%timeit expr = oe.contract_expression(einsum_str, *shapes, memory_limit=-1)
#> 76.2 ms ± 1.05 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)


This is pretty manageable, though we might want to think about splitting the expression up if we go a lot bigger. Importantly, we can then use this repeatedly with any set of matching arrays:

arrays = [np.random.randn(*shp) / 4 for shp in shapes]
expr(*arrays)
#> array(23.23628116)

arrays = [np.random.randn(*shp) / 4 for shp in shapes]
expr(*arrays)
#> array(-12.21091879)


### Full path

And if we really want we can generate the full contraction path info:

print(oe.contract_path(einsum_str, *arrays, memory_limit=-1)[1])
#>          Naive scaling:  298
#>      Optimized scaling:  5
#>       Naive FLOP count:  1.031e+248
#>   Optimized FLOP count:  1.168e+06
#>    Theoretical speedup:  88264689284468460017580864156865782413140936705854966013600065426858041248009637246968036807489558012989638169986640870276510490846199301907401763236976204166215471281505344088317454144870323271826022036197984172898402324699098341524952317952.000
#>   Largest intermediate:  3.000e+02 elements
#> --------------------------------------------------------------------------------
#> scaling        BLAS                current                             remaining
#> --------------------------------------------------------------------------------

Where we can see the speedup over a naive einsum is about 10^241, not bad!