# 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 get_symbol().

Let’s set up the required einsum 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)


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

...

4           TDOT            Ğğ,ĠğĢ->ĠĞĢ                  ĠĞġ,ģĢĥ,ģġĤ,Ĥĥ,ĠĞĢ->
4           GEMM            ĠĞĢ,ĠĞġ->ġĢ                       ģĢĥ,ģġĤ,Ĥĥ,ġĢ->
4           GEMM            Ĥĥ,ģĢĥ->ģĢĤ                          ģġĤ,ġĢ,ģĢĤ->
4           TDOT            ģĢĤ,ģġĤ->ġĢ                               ġĢ,ġĢ->
2            DOT                ġĢ,ġĢ->                                    ->


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