Repository: dgasmith/opt_einsum
Pull request #102 (2019-08-28 11:45:00)

added dynamic programming path
## Description
In the field of tensor network states (matrix product states, projected entangled pair states, ...) the most common approach to optimise tensor network contractions is based on dynamic programming [Phys. Rev. E 90, 033315 (2014)], which usually provides an optimal contraction path within a reasonable amount of time (even for a large number of tensors). This pull request provides such an implementation.

## Questions
- [ ] Is there a simple way to run the benchmarks of the "Performance Comparison" section of the documentation for the dynamic programming optimiser to compare with the other approaches?

codecov[bot]: # [Codecov](https://codecov.io/gh/dgasmith/opt_einsum/pull/102?src=pr&el=h1) Report
> :exclamation: No coverage uploaded for pull request base (`master@f0ec70f`). [Click here to learn what that means](https://docs.codecov.io/docs/error-reference#section-missing-base-commit).
> The diff coverage is `14.28%`.


codecov[bot]: # [Codecov](https://codecov.io/gh/dgasmith/opt_einsum/pull/102?src=pr&el=h1) Report
> :exclamation: No coverage uploaded for pull request base (`master@f0ec70f`). [Click here to learn what that means](https://docs.codecov.io/docs/error-reference#section-missing-base-commit).
> The diff coverage is `14.28%`.


jcmgray: > Do you have standard test cases you are using for the other optimisers as well to ensure they are yielding correct results?

Yes, check the ``test_contract.py`` and ``test_paths.py`` files - should just be a matter of adding a few parametrizations to the existing tests.

Those tweaks seem to have sped it up a fair bit - if you have any other potential optimizations in mind let me know, else I will run the benchmarks in a bit!
mrader1248: This code should work for contractions with summation indices occurring three or more times (although I am not entirely sure, how optimal the resulting contractions will be in general) and it also works for non-summation indices occuring more than once.

Do you have standard test cases you are using for the other optimisers as well to ensure they are yielding correct results?

It would be great, if you could run the performance comparison.
jcmgray: Not quite sure what you mean by 'repair', ``"ab,ab,c"`` would generate the output automatically as ``"ab,ab,c->c"``, but when the arrow is explicitly there ``"ab,ab,c->"`` it implies that the ``'c'`` dimension is indeed summed over and won't appear in the output, so this is still a valid contraction.
mrader1248: I added `'dp'` to `all_optimizers` to test my implementation. You can still wait with running the benchmark until there are no more failures.
Something I just realised: The code in `contract` does not 'repair' the output indices before calling an optimiser? At the moment test cases like `"ab,ab,c->"` fail for my implementation as it assumes complete specification of the output indices (it works for `"ab,ab,c->c"` and `"ab,ab,c->abc"`). Do you have some common routine to repair the output indices?
jcmgray: Here's the initial results I have for path finding time:

![oe-dp-times](https://user-images.githubusercontent.com/8982598/63879441-a543cb00-c9c3-11e9-9897-7d09ebf90b2c.png)

And FLOPS improvement over ``greedy``:

![oe-dp-flops](https://user-images.githubusercontent.com/8982598/63879462-af65c980-c9c3-11e9-9163-fb92aa3deaf1.png)

For this class of problems (basically random regular like contractions) it seems to give similar quality to ``branch-all`` but is faster for 10+ tensors.
mrader1248: Thank you for explaining this example to me. Indeed, I thought `"ab,ab,c"` and `"ab,ab,c->"` is the same. With 'repair' I meant that `contract` should transform `"ab,ab,c"` into `"ab,ab,c->c"` before passing it to the optimiser (and I just realised that this is actually done). So, nevermind.

Indeed, my code fails for summation indices occurring only once. However, as far as I understand it, summing over these indices in the beginning is always optimal?
mrader1248: > When you say "summation indices occurring only once" do you refer to `"ab,ab,c"` like operations in all cases? If so, I believe this is true and something that we have been toying around with as well. See [here](https://github.com/dgasmith/opt_einsum/issues/99).
> 
> Not exactly identical, but similar problems. I do believe completing all Hadamard products and forming all scalar values first will always be an optimal strategy.

It seems very reasonable that first summing over indices occurring only in one input (e.g. `"c"` and `"d"` for `"ab,ab,c,dd->"`) is indeed optimal (at least I cannot think of any counter example). I put the code achieving this in a function `_sum_single_tensor_indices(inputs, output)` if anybody wants to use it.

@jcmgray: Thank you for these first benchmarks. I realised there might be one more optimisation I should implement (decomposing the contraction graph into disconntected subgraphs first).
mrader1248: I added another improvement, which causes the optimiser to scale only linearly with the number of disconnected subgraphs of the contraction graph. If all the tests are successful for the latest commit, it would be interesting to see if/how performance improved.

So far I only added the `'dp'` path to the tests in `test_contract.py`. What kind of tests have to be added in `test_paths.py`?
mrader1248: > With regard to `test_paths.py`, maybe just add it to the `path_edge_tests` cases to make sure it matches `optimal` and maybe to `test_can_optimize_outer_products` as well.

The problem with `path_edge_tests` is that the solutions are not unique. Considering the first test `'eb,cb,fb->cef'` (apparently it's tested with shapes `[(4, 3), (4, 3), (3, 3)]`) it doesn't matter if the contraction sequence is `((0, 2), (0, 1))` or `((1, 2), (0, 1))` (all your optimisers yield the first solution, my optimiser yields the second one).

The test `test_can_optimize_outer_products` has the same issue: The result is not unique. For `"ab,cd,ef,fg"` it doesn't matter if the summation over `f` is performed before or after the first outer product. However, at this point I realised something different as well: For an outer product of three tensors, I thought the sequence `[(0,1,2)]` should be used instead of `[(0,1), (0,1)]`. However, the latter version is twice as fast. Do you know why?

> Might it also be possible to lint it for PEP8 as well? (@dgasmith not sure if this is official standard for `opt_einsum` but certainly the rest of the codebase seems to conform currently).

Can you tell me where I violate PEP8?
jcmgray: > The problem with path_edge_tests is that the solutions are not unique. Considering the first test 'eb,cb,fb->cef' (apparently it's tested with shapes [(4, 3), (4, 3), (3, 3)]) it doesn't matter if the contraction sequence is ((0, 2), (0, 1)) or ((1, 2), (0, 1)) (all your optimisers yield the first solution, my optimiser yields the second one).

Yes I think these can just be thought of as regression tests, I'd just put in the path that ``'dp'`` currently finds, but personally also happy if you want to skip. By the way, is the path finding deterministic?

> For an outer product of three tensors, I thought the sequence [(0,1,2)] should be used instead of [(0,1), (0,1)]. However, the latter version is twice as fast. Do you know why?

I think in general, a pairwise contraction sequence always is better. For the first path you have ``2 * d^3`` operations, and in the second case ``d^2 + d^3`` operations. Moreover, in practical cases, the pairwise paths can use tensordot and BLAS, whereas if there are more than 2 tensors, ``einsum`` must be used.

> Can you tell me where I violate PEP8?

There are quite a few, but potentially just ones that are turned off in your editor! Here's what I get:
```
paths.py:16:73: W291 trailing whitespace
paths.py:662:6: W291 trailing whitespace
paths.py:669:6: W291 trailing whitespace
paths.py:670:79: W291 trailing whitespace
paths.py:673:1: W293 blank line contains whitespace
paths.py:676:1: W293 blank line contains whitespace
paths.py:677:12: E261 at least two spaces before inline comment
paths.py:680:1: W293 blank line contains whitespace
paths.py:684:1: W293 blank line contains whitespace
paths.py:688:1: W293 blank line contains whitespace
paths.py:692:1: W293 blank line contains whitespace
paths.py:697:1: W293 blank line contains whitespace
paths.py:700:1: W293 blank line contains whitespace
paths.py:703:1: W293 blank line contains whitespace
paths.py:706:20: E272 multiple spaces before keyword
paths.py:706:50: E271 multiple spaces after keyword
paths.py:708:1: W293 blank line contains whitespace
paths.py:717:1: W293 blank line contains whitespace
paths.py:725:1: W293 blank line contains whitespace
paths.py:726:60: E226 missing whitespace around arithmetic operator
paths.py:728:1: W293 blank line contains whitespace
paths.py:729:40: E261 at least two spaces before inline comment
paths.py:730:1: W293 blank line contains whitespace
paths.py:741:1: W293 blank line contains whitespace
paths.py:743:1: W293 blank line contains whitespace
paths.py:748:1: W293 blank line contains whitespace
paths.py:752:1: W293 blank line contains whitespace
paths.py:755:1: W293 blank line contains whitespace
paths.py:759:22: E261 at least two spaces before inline comment
paths.py:762:1: W293 blank line contains whitespace
paths.py:765:1: W293 blank line contains whitespace
paths.py:771:1: W293 blank line contains whitespace
paths.py:773:76: W291 trailing whitespace
paths.py:774:1: W293 blank line contains whitespace
paths.py:775:116: W291 trailing whitespace
paths.py:783:54: E261 at least two spaces before inline comment
paths.py:785:81: E261 at least two spaces before inline comment
paths.py:790:1: W293 blank line contains whitespace
paths.py:792:1: W293 blank line contains whitespace
paths.py:802:79: W291 trailing whitespace
paths.py:803:63: W291 trailing whitespace
paths.py:804:77: W291 trailing whitespace
paths.py:806:71: W291 trailing whitespace
paths.py:808:1: W293 blank line contains whitespace
paths.py:810:71: W291 trailing whitespace
paths.py:811:72: W291 trailing whitespace
paths.py:849:1: W293 blank line contains whitespace
```
mrader1248: > Yes I think these can just be thought of as regression tests, I'd just put in the path that `'dp'` currently finds, but personally also happy if you want to skip. By the way, is the path finding deterministic?

Yes, the path finding is deterministic. I personally just don't like these kind of tests: If somebody changes the path finding function in the future in a way that will affect the resulting path, it will take an annoying amount of time to find out that actually the test requires a specific solution.

> I think in general, a pairwise contraction sequence always is better. For the first path you have `2 * d^3` operations, and in the second case `d^2 + d^3` operations. Moreover, in practical cases, the pairwise paths can use tensordot and BLAS, whereas if there are more than 2 tensors, `einsum` must be used.

Why is it `2 * d^3`? I thought, it should be `d^3` (looping over three indices and multiplying three numbers each iteration)? However, I see that for pairwise outer products you can use BLAS and BLAS of course always wins. ;-)

Ok, I will take a look at these PEP8 violations.
jcmgray: > Yes, the path finding is deterministic. I personally just don't like these kind of tests: If somebody changes the path finding function in the future in a way that will affect the resulting path, it will take an annoying amount of time to find out that actually the test requires a specific solution.

Yes I can appreciate this reasoning - I guess ideally we would test against equivalent paths, but for the moment best to leave.

> Why is it 2 * d^3? I thought, it should be d^3 (looping over three indices and multiplying three numbers each iteration)? However, I see that for pairwise outer products you can use BLAS and BLAS of course always wins. ;-)

There is some subtlety here to do with compiler instructions etc but essentially ``opt_einsum`` assumes that multiplying 3 numbers is twice as expensive as multiplying two numbers (i.e. ``x * y * z == (x * y) * z``). Equivalently, if done that way ``(x * y)`` is recomputed for each ``d`` times redundantly.
mrader1248: I fixed the bug I mentioned in the context of `_sum_single_tensor_indices`. This also led to the removal of this function.

I added the dp path to the tests in `test_paths.py` and used the results, my optimiser currently finds (and checked by hand that they are correct).

The performance of my optimiser might have slightly improved as I removed two set intersections from a condition as they are no longer needed due to the initial decomposition into disonnected subgraphs.

`pycodestyle --max-line-length=120 paths.py` now gives me no more errors.

One more thing I found out: Running the rather trivial contraction code
```
expr = oe.contract_expression("i,i,j,j,k,l->", *[v.shape]*6, optimize=...)
%timeit expr(*[v]*6)
```
gives me a runtime of ~2.2 ms for the optimal and ~1.3 ms for the dp path (for v with 1M float64 elements) and ~9 s (!) for the greedy path (for only 100k elements). The optimal path is `[(4,), (4,), (4, 5), (0, 1), (2, 3), (0, 1), (0, 1)]` and my dp path is `[(4,), (4,), (4, 5), (0, 1), (2, 3), (0, 1), (0, 1)]`. I don't really understand how my code can be faster than the optimal path (especially because the dp path has 7 instead of 5 elements and both should have require `4d` operations), but I thought this might be interesting for you.
Although this is quite off-topic now: The greedy path produces a `d^2`-solution which could be easily circumvented by performing all the summations over single inputs in the very beginning.
jcmgray: Amazing, thanks for those changes. I will try and re-run the benchmarks soon, probably with all the other optimizers actually again since I have changed machines.

> I don't really understand how my code can be faster than the optimal path (especially because the dp path has 7 instead of 5 elements and both should have require 4d operations), but I thought this might be interesting for you.

The estimated flop cost, while well correlated with actual run time, is still just an estimate. There are many factors to do with copying, transposing, threadability ... etc that means some paths are just noticeably faster despite the same estimated cost.

> Although this is quite off-topic now: The greedy path produces a d^2-solution which could be easily circumvented by performing all the summations over single inputs in the very beginning.

I'm not sure I follow this, which bit of greedy? As far as I am aware the single d^2 step of greedy is to find all the initial inner products.
mrader1248: > Amazing, thanks for those changes. I will try and re-run the benchmarks soon, probably with all the other optimizers actually again since I have changed machines.

Great. I am looking forward to the benchmark results.

> I'm not sure I follow this, which bit of greedy? As far as I am aware the single d^2 step of greedy is to find all the initial inner products.

With `d` I meant the length of `v`. The greedy path is `[(0, 1), (0, 1), (0, 1), (0, 1), (0, 1)]`. The first two operations are the two scalar products [`O(d)`] and the third operation is a sum over an outer product [`O(d^2)`].
jcmgray: > With d I meant the length of v. The greedy path is [(0, 1), (0, 1), (0, 1), (0, 1), (0, 1)]. The first two operations are the two scalar products [O(d)] and the third operation is a sum over an outer product [O(d^2)].

Ah sorry yes I see. At least personally these trivial cases of disconnected contractions haven't been on my radar which is maybe also why they aren't handled sensibly in various places. Combined with #99 it might be time to introduce some universal preprocessing of all the trivial contractions like these + hadamard products etc.

---

Here are the results, which look good! I think the apparent improvement over previous reuslts is mostly due to running with the full set of contractions as the others which I had not checked (sorry):

![oe-dp-times](https://user-images.githubusercontent.com/8982598/64038182-5a56be80-cb4f-11e9-8376-10f4cb4b53cb.png)

![oe-dp-flops](https://user-images.githubusercontent.com/8982598/64038190-5d51af00-cb4f-11e9-8013-4f880ce066bf.png)

Certainly this method seems to outperform ``branch-all``. I've also included ``random-greedy-128`` for reference as I think it might be nice to have an ``'auto-hq'`` mode that targets a longer solving time first using ``'dp'`` before switching over to maybe ``'random-greedy'``.

---

Just out of interest, can the ``dp`` approach be truncated and so made anytime? I.e. is there way to sort the search such that likely contractions are checked first or does it have to be exhaustive? Otherwise I'm thinking the ``BranchBound`` path finders might be refitted to fill that niche. In some brief tests ``branch-all`` performs very nicely if you just make it timeout rather than checking every solution. [No worries if you don't want to think about these things :)].
mrader1248: Thank you very much for rerunning the benchmarks. It's very nice to see that `dp` is speedup-wise on par with the `optimal` approach, but can be used for even larger contractions.

Unfortunately, I don't see a trivial way to run a truncated version of `dp`. The `dp` optimiser simply finds the optimal contractions for all sets of n (connected) tensors by combining optimal contractions of all disjoint sets of m and n - m (connected) tensors. So I'm not sure if any other optimiser could use these intermediate results.

I implemented two more small code optimisations that should make `dp` run ~20-25 % faster. Before doing anything, my code converts all the n occurring indices to integers 0, 1, ..., n-1. I guess this could also be used to speed up other optimisers (e.g. operations on sets of integers are faster in CPython than on sets of strings). Another optimisation was to use bitmaps (integers) instead of tensor sets. Furthermore, FYI: I observed that ~40 % of `dp` runtime is spent in `compute_size_by_dict`.
dgasmith: Not sure I have seen a way that speed up `dp` algorithms either.

The bitwise operations are interesting, since the `dp` algorithm can only handle pretty tame numbers of total indices this is likely fine for most cases. We might want to think if there are bounds in the number of supported input tensors and indices.

`compute_size_by_dict` is a good candidate to cythonize first. If we remap the code base to store sizes as an array and indices become integers that map that array we could make several items more efficient. 
mrader1248: As far as I understand it, this approach should be correct for an arbitrarily large number of inputs, as python automatically switches to some kind of big integer type. However, when using this big integer type, it could slow down these bitwise operations (which might still be faster than using the set classes).

Yes, to cythonize `compute_size_by_dict` might be a good idea, although this is beyond the scope of this pull request. However, keep in mind that I'm now calling this function with an array (=dict mapping integer indices to sizes) instead of a dict.
jcmgray: Another update to the benchmarks then! By the way I've updated the notebook [here](https://gist.github.com/jcmgray/702a70c188b20e74506e635da2033ab0) that can generate these plots.

![oe-times-graph](https://user-images.githubusercontent.com/8982598/64082414-7baddb00-cd06-11e9-8863-a5a5b6cec9cc.png)

![oe-flops-graph](https://user-images.githubusercontent.com/8982598/64082415-7e103500-cd06-11e9-84a4-5c38edfc2e6e.png)

Have you got any other planned updates @mrader1248? I'm thinking ``'dp'`` should probably replace ``'branch-all'`` in the ``'auto'`` path for ``len(inputs)=[5, 6]``, and at some point it would be good to add the dynamic programming approach to the proper documents, but maybe these could be done in a later, separate PR if you are keen to get this in?
mrader1248: No, I'm fine now with my code. So I'm not planning to change anything else at the moment. I feel honoured that you want to put `dp` immediately into the `auto`-chain, but I'm not sure whether this should be accomplished with my pull request. The choice which optimisers should be part of the `auto`-chain is completely yours.

What to you mean with 'the proper documents'?
jcmgray: > What to you mean with 'the proper documents'?

``opt_einsum`` has documentation e.g. [on the path finders](https://optimized-einsum.readthedocs.io/en/latest/path_finding.html) that is built from the files in ``opt_einsum/docs/source``, see for example [the page on the branching path](https://github.com/dgasmith/opt_einsum/blob/master/docs/source/branching_path.rst). At some point it would be good to include information on this new path.

> The choice which optimisers should be part of the auto-chain is completely yours.

@dgasmith what do you think? It seems like a fairly natural upgrade to me.
dgasmith: I agree on replacing the branching algorithm here with ‘dp’, the statistics look good. At some point we are going to have to do a white paper on these algorithms. We are collecting enough it could be interesting. 
mrader1248: I remembered that a good upper bound for the contraction cost could give quite some speedup, so I added an optional parameter `cost_limit` for `dp` to reject contractions (already of subsets of the tensors) that exceed this limit. Then I also realised that for a larger number of inputs, I could use an estimate from `random-greedy-128` to speed up `dp`.
jcmgray: Sounds like a good enhancement. Rather than adding a new cost function how possible would it be to make the cost the same (if it is not currently) as the flops estimate as used by the other contractors? It would be also nice to add ``cost_limit`` to the constructor of ``DynamicProgrammingOptimizer``. That would allow this kind of thing:

```python
rg_opt = oe.RandomGreedy()
oe.contract_path(eq, *arrays, optimize=rg_opt)

dp_opt = oe.DynamicProgrammingOptimizer(cost_limit=rg_opt.best['flops'])
oe.contract_path(eq, *arrays, optimize=dp_opt)
```
And does running just a simple `greedy` cost_limit for all sizes help?
mrader1248: I'm not sure, how this flops estimate is determined, but it seems to be larger than the one I'm using (I want to have it as small as possible, but still large enough).

You are right, I moved `cost_limit` to the constructor.

I guess this should help quite a bit. Just take this example of 20 inputs (5x4 tensors connected on a square lattice):
```
inputs = [
    {0, 15}, {0, 1, 16}, {1, 2, 17}, {2, 18},
    {3, 15, 19}, {3, 4, 16, 20}, {4, 5, 17, 21}, {5, 18, 22},
    {6, 19, 23}, {6, 7, 20, 24}, {7, 8, 21, 25}, {8, 22, 26},
    {9, 23, 27}, {9, 10, 24, 28}, {10, 11, 25, 29}, {11, 26, 30},
    {12, 27}, {12, 13, 28}, {13, 14, 29}, {14, 30}
]
size_dict = {k: 2 for k in set.union(*inputs)}
```
`random-greedy-128` takes 157 ms ± 7.12 ms, `dp` without `cost_limit` ~ 100 s and `dp` with the `random-greedy-128` estimate required 563 ms ± 2.44 ms in my test.
jcmgray: > I'm not sure, how this flops estimate is determined, but it seems to be larger than the one I'm using (I want to have it as small as possible, but still large enough).

Yes I mean ideally ``'dp'`` would use this same flop costing internally as well, based on calling ``helpers.flop_count``, then *any* cost calculated by another path would be a valid sieve for this.
mrader1248: Oh, now I see how you meant it. Of course, this I could do. However, I don't quite understand the `flop_count` function. First of all, executing the examples from the docstring yield other results than stated in the docstring. Secondly, what is this `inner` parameter for?
jcmgray: Yes I was confused a bit about [``flop_count`` and inner here](https://github.com/dgasmith/opt_einsum/pull/73#issuecomment-434897041), but as @dgasmith pointed out, it's because if there is an inner product, each element involves both a multiplication *and* a summation (subtlety ... unless the cpu instructions are executed with FMA support).

Regarding the docstring examples, looks like they are out-of-date and should be corrected.

At some point a possibility would be to have the current advanced flop estimate changeable so as to allow for instance the cost you are currently using instead. I think these are essentially the same (often just differ by a factor of 2) as ``'dp'`` currently finds the same paths as ``'optimal'``.
mrader1248: If I just contract two tensors, `inner` has to be true if it is a scalar product, i.e. if there are no more indices left after performing this contraction? In the `dp` code I'm using now:
```
cost = cost1 + cost2 + helpers.flop_count(i1_union_i2, len(i) == 0, 2, size_dict)
```
The ugly thing about it: I have to compute `i` now before estimating the cost, which makes it more expensive..

But what's anyway quite strange: To get an upper bound, I use
```
rg = path_random.RandomGreedy(max_repeats=128)
p = rg(inputs, output, size_dict, memory_limit)
cost_limit = 1.1 * 0.5 * rg.best["flops"]
```
I can halve the flops from `random-greedy-128` and still get all the tests correct..

Also, I redefined the function `path_cost` I introduced for cross checking:
```
def path_cost(inputs, output, size_dict, path):
    c = 0
    for p in path:
        inputs_to_contract = [i for j, i in enumerate(inputs) if j in p]
        i_inputs_to_contract = set.union(*inputs_to_contract)
        inputs = [i for j, i in enumerate(inputs) if j not in p]
        i_contract = i_inputs_to_contract - output
        if len(inputs) > 0:
            i_contract -= set.union(*inputs)
        new_input = i_inputs_to_contract - i_contract
        c += flop_count(i_inputs_to_contract, len(new_input) == 0, len(inputs_to_contract), size_dict)
        inputs += [new_input]
    return c
```
If I use this function to recompute the cost of the `random-greedy-128` path, it's only half as much as `rg.best["flops"]`. So, is this function wrong or is there a bug in the flops estimate of `random-greedy-128`?
jcmgray: Yes I think you need to tweak which indices are supplied. Here's a version that works for flops:

```python
def path_cost_flops(inputs, output, size_dict, path):
    inputs = list(inputs)
    c = 0
    for p in path:
        inputs_to_contract = [inputs.pop(i) for i in sorted(p, reverse=True)]
    
        either = set.union(*inputs_to_contract)
        shared = set.intersection(*inputs_to_contract)
        keep = set.union(output, *inputs)
    
        c += helpers.flop_count(either, shared - keep, len(p), size_dict)
        inputs.append(either & keep)
        
    return c
```
If it is significantly more expensive to keep track of flops like this then maybe sticking to the original cost is the best thing to do. I might actually run some timings to see how well each estimated cost correlates with real contraction time.
jcmgray: Getting slightly off track here but this is quite interesting. If you generate random expressions, and compute the full flop cost (as currently done) and the 'number of operations' cost, then compared to the actual time to perform the computation, the slightly simpler 'number of operations' does seem like a slightly better predictor:

First is histogram of ``log10(flop_cost / time)``.
Second is histogram of ``log10(2 * num_ops / time)``.
Both are coloured by number of output indices.

![visualization(4)](https://user-images.githubusercontent.com/8982598/64127211-9df83500-cda8-11e9-9db8-9eb7bc237269.png)

Both predict the actual time will be about 1e-8 times the cost, and they match when ``n_out == 0``. The second is very slightly tighter but also has no dependency on ``n_out``.

All of which is to say, it would not necessarily be a bad thing to simplify the FLOP counting across the board.
mrader1248: I do not have high precision measurements, but `dp` is certainly 30-40 % slower when I have to compute the remaining indices of a contraction beforehand. So from this point of view it might really make sense to stick to the current estimators.

Could you be so kind as to rerun the `dp` benchmarks? In contrast to the previous changes I guess these latest changes should really affect the optimisation times.
jcmgray: I'll finish running the updated benchmarks soon - as @mrader1248 says, implementing the cost limit with a greedy call should make a big difference! But otherwise agree with leaving the other stuff to different issues.
mrader1248: I made a last modification: Instead of running random-greedy in the beginning, I am now using an iteratively increasing cost cap. For the example of 5x4 tensors on a square lattice this brings me down from 563 ms to 420 ms and I can do 6x5 tensors in less than 30 s. The problem with the cost estimate from random-greedy is that for a larger number of inputs it becomes too bad to really speed up dp. Of course, also with the iterative cost cap dp scales exponentially with the number of inputs, but the overall performance is on the one hand side a bit better and on the other hand side the dependency on random-greedy is removed.

@dgasmith, I also think -- if you are okay with the latest version of my pull request -- that it is time to merge. Regarding the points you mentioned above: I think that these four points should be tackled in new issues/pull requests.
jcmgray: Happy to merge if you want to get this in - not trying to delay! 

However I think the cost cap stuff does merit a bit more thinking about maybe in a later PR as there likely gains to be had. For any planar graphs (such as lattices) the random-greedy approach generally finds the optimal path, indeed, I find for both 4x5 and 5x6 lattices it finds the exact same path cost as ``'dp'``. So it seems surprising it wasn't working as a good cost cap. Out of interest, why did you opt for ``int(1.1 * cost_limit)``, surely any valid path cost found is a valid limit?

Here is a graph of timings including the newest method (``'dp-costcap``), and I ran the greedy labelled ones with that 1.1 scaling removed:

![oe-cost-time](https://user-images.githubusercontent.com/8982598/64253131-7faf4800-cf14-11e9-8e88-4a37cbd66734.png)

It might be worth in the end having options for:
* using no costlimit
* using the iterative method
* supplying a str or PathOptimizer to precompute one
* using some automatically selected number of greedy runs
dgasmith: Tiers of greedy runs could be the name of the day:
```python
ngreedy = 0
if nterms > 7:
   ngreedy = (terms - 7) ** 2
ngreedy = max(ngreedy, 128)
```

Though simply `greedy-1` seems to do a reasonable job across the board.

Not trying to push this through too quickly. Just wanted to rule a few items that are in and out of scope of this PR so that it doesn't go on forever :)
mrader1248: Maybe it's really the best to just merge this for the moment to have dp in opt_einsum and then start  thinking about possibilities of using combinations of dp and random-greedy. I would go for the version with the cost cap as it's completely self-contained and the performance is not bad.

@jcmgray, I used the factor `1.1` as I wasn't entirely sure whether there could be a small +epsilon in the cost. But of course, if this isn't the case one can omit this factor.

@dgasmith, I already realised that with an increasing number of inputs also random-greedy-X with larger X should be used. However, I have the bad feeling that the optimal X might also scale exponentially with the number of inputs.

Do you want to have any further changes in this PR or are you okay with merging it?
dgasmith: Well, not sure if we are really after *optimal* X, but just something reasonable at first.

Looking over again, now. But overall LGTM.

@jcmgray Are you ok with merging?
jcmgray: Yes happy to merge as is - I've updated that latest graph with ``n=20`` and the various methods look pretty close in the limit. Presumably those last few review comments should be addressed first, then yes all looks good to me!
mrader1248: > Well, not sure if we are really after _optimal_ X, but just something reasonable at first.

The big advantage of the cost capping strategy is that we don't need to worry about this X.

@dgasmith, I incorporated your suggestions regarding documentation. From my point of view we are now ready to merge.

Patch:
diff --git a/opt_einsum/paths.py b/opt_einsum/paths.py
index 0a18b65..697f866 100644
--- a/opt_einsum/paths.py
+++ b/opt_einsum/paths.py
@@ -12,7 +12,10 @@ import numpy as np
 
 from . import helpers
 
-__all__ = ["optimal", "BranchBound", "branch", "greedy", "auto", "get_path_fn"]
+__all__ = [
+    "optimal", "BranchBound", "branch", "greedy", "auto", "get_path_fn",
+    "DynamicProgrammingOptimizer", "dynamic_programming"
+]
 
 
 _UNLIMITED_MEM = {-1, None, float('inf')}
@@ -652,6 +655,324 @@ def greedy(inputs, output, size_dict, memory_limit=None, choose_fn=None, cost_fn
     return ssa_to_linear(ssa_path)
 
 
+def _tree_to_sequence(c):
+    """
+    Converts a contraction tree to a contraction path as it has to be
+    returned by path optimizers. A contraction tree can either be an int
+    (=no contraction) or a tuple containing the terms to be contracted. An
+    arbitrary number (>= 1) of terms can be contracted at once. Note that
+    contractions are commutative, e.g. (j, k, l) = (k, l, j). Note that in
+    general, solutions are not unique.
+
+    Parameters
+    ----------
+    c : tuple or int
+        Contraction tree
+
+    Returns
+    -------
+    path : list[set[int]]
+        Contraction path
+
+    Examples
+    --------
+    >>> _tree_to_sequence(((1,2),(0,(4,5,3))))
+    [(1, 2), (1, 2, 3), (0, 2), (0, 1)]
+    """
+
+    # ((1,2),(0,(4,5,3))) --> [(1, 2), (1, 2, 3), (0, 2), (0, 1)]
+    #
+    # 0     0         0           (1,2)       --> ((1,2),(0,(3,4,5)))
+    # 1     3         (1,2)   --> (0,(3,4,5))
+    # 2 --> 4     --> (3,4,5)
+    # 3     5
+    # 4     (1,2)
+    # 5
+    #
+    # this function iterates through the table shown above from right to left;
+
+    if type(c) == int:
+        return []
+
+    c = [c]  # list of remaining contractions (lower part of columns shown above)
+    t = []   # list of elementary tensors (upper part of colums)
+    s = []   # resulting contraction sequence
+
+    while len(c) > 0:
+        j = c.pop(-1)
+        s.insert(0, tuple())
+
+        for i in sorted([i for i in j if type(i) == int]):
+            s[0] += (sum(1 for q in t if q < i),)
+            t.insert(s[0][-1], i)
+
+        for i in [i for i in j if type(i) != int]:
+            s[0] += (len(t) + len(c),)
+            c.append(i)
+
+    return s
+
+
+def _find_disconnected_subgraphs(inputs, output):
+    """
+    Finds disconnected subgraphs in the given list of inputs. Inputs are
+    connected if they share summation indices. Note: Disconnected subgraphs
+    can be contracted independently before forming outer products.
+
+    Parameters
+    ----------
+    inputs : list[set]
+        List of sets that represent the lhs side of the einsum subscript
+    output : set
+        Set that represents the rhs side of the overall einsum subscript
+
+    Returns
+    -------
+    subgraphs : list[set[int]]
+        List containing sets of indices for each subgraph
+
+    Examples
+    --------
+    >>> _find_disconnected_subgraphs([set("ab"), set("c"), set("ad")], set("bd"))
+    [{0, 2}, {1}]
+
+    >>> _find_disconnected_subgraphs([set("ab"), set("c"), set("ad")], set("abd"))
+    [{0}, {1}, {2}]
+    """
+
+    subgraphs = []
+    unused_inputs = set(range(len(inputs)))
+
+    i_sum = set.union(*inputs) - output  # all summation indices
+
+    while len(unused_inputs) > 0:
+        g = set()
+        q = [unused_inputs.pop()]
+        while len(q) > 0:
+            j = q.pop()
+            g.add(j)
+            i_tmp = i_sum & inputs[j]
+            n = {k for k in unused_inputs if len(i_tmp & inputs[k]) > 0}
+            q.extend(n)
+            unused_inputs.difference_update(n)
+
+        subgraphs.append(g)
+
+    return subgraphs
+
+
+def _bitmapset_indices(s):
+    """
+    Returns a generator object allowing to iterate over the elements contained
+    in a bitmap set.
+
+    Parameters
+    ----------
+    s : int
+        The bitmap set to iterate over
+
+    Returns
+    -------
+    path : generator
+        Generator object to iterate over the elements in s
+
+    Examples
+    --------
+    >>> type(_bitmapset_indices(0b1001011))
+    generator
+
+    >>> list(_bitmapset_indices(0b1001011))
+    [0, 1, 3, 6]
+    """
+    j = 0
+    while s != 0:
+        if s & 1 != 0:
+            yield j
+        s >>= 1
+        j += 1
+
+
+class DynamicProgrammingOptimizer(PathOptimizer):
+    """
+    Finds the optimal path of pairwise contractions without intermediate outer
+    products based a dynamic programming approach presented in
+    Phys. Rev. E 90, 033315 (2014) (the corresponding preprint is publically
+    available at https://arxiv.org/abs/1304.6112). This method is especially
+    well-suited in the area of tensor network states, where it usually
+    outperforms all the other optimization strategies.
+
+    This algorithm shows exponential scaling with the number of inputs
+    in the worst case scenario (see example below). If the graph to be
+    contracted consists of disconnected subgraphs, the algorithm scales
+    linearly in the number of disconnected subgraphs and only exponentially
+    with the number of inputs per subgraph.
+    """
+
+    def __call__(self, inputs, output, size_dict, memory_limit=None):
+        """
+        Parameters
+        ----------
+        inputs : list
+            List of sets that represent the lhs side of the einsum subscript
+        output : set
+            Set that represents the rhs side of the overall einsum subscript
+        size_dict : dictionary
+            Dictionary of index sizes
+        memory_limit : int
+            The maximum number of elements in a temporary array
+
+        Returns
+        -------
+        path : list
+            The contraction order (a list of tuples of ints).
+
+        Examples
+        --------
+        >>> n_in = 3  # exponential scaling
+        >>> n_out = 2 # linear scaling
+        >>> s = dict()
+        >>> i_all = []
+        >>> for _ in range(n_out):
+        >>>     i = [set() for _ in range(n_in)]
+        >>>     for j in range(n_in):
+        >>>         for k in range(j+1, n_in):
+        >>>             c = oe.get_symbol(len(s))
+        >>>             i[j].add(c)
+        >>>             i[k].add(c)
+        >>>             s[c] = 2
+        >>>     i_all.extend(i)
+        >>> o = DynamicProgrammingOptimizer()
+        >>> o(i_all, set(), s)
+        [(1, 2), (0, 4), (1, 2), (0, 2), (0, 1)]
+        """
+
+        # convert all indices to integers (makes set operations ~10 % faster)
+        symbol2int = {c: j for j, c in enumerate(set.union(*inputs) | output)}
+        inputs = [set(symbol2int[c] for c in i) for i in inputs]
+        output = set(symbol2int[c] for c in output)
+        size_dict = {symbol2int[c]: v for c, v in size_dict.items() if c in symbol2int}
+        size_dict = [size_dict[j] for j in range(len(size_dict))]
+
+        # all summation indices occurring exactly in one input:
+        i_single = set(
+            c for c in set.union(*inputs) - output
+            if sum(1 for i in inputs if c in i) == 1
+        )
+
+        # contraction expressions for all inputs that have already been
+        # reduced to scalars:
+        inputs_done = [
+            (j,) for j, i in enumerate(inputs)
+            if len(i - i_single) == 0
+        ]
+
+        # remaining input index sets and corresponding contraction expressions;
+        # indices from i_single are removed and if a single-tensor contraction
+        # is performed, the contraction expression is (j,) instead of j;
+        inputs, inputs_contractions = zip(*[
+            (i - i_single, j if i.isdisjoint(i_single) else (j,))
+            for j, i in enumerate(inputs)
+            if len(i - i_single) > 0
+        ])
+
+        # a list of all neccessary contraction expressions for each of the
+        # disconnected subgraphs and their size
+        subgraph_contractions = inputs_done
+        subgraph_contractions_size = [1]*len(inputs_done)
+
+        for g in _find_disconnected_subgraphs(inputs, output):
+
+            # dynamic programming approach to compute x[n] for subgraph g;
+            # x[n][set of n tensors] = (indices, cost, contraction)
+            # the set of n tensors is represented by a bitmap: if bit j is 1,
+            # tensor j is in the set, e.g. 0b100101 = {0,2,5}; set unions
+            # (intersections) can then be computed by bitwise or (and);
+            x = [None]*2 + [dict() for j in range(len(g)-1)]
+            x[1] = {1 << j: (inputs[j], 0, inputs_contractions[j]) for j in g}
+
+            # convert set of tensors g to a bitmap set:
+            g = functools.reduce(lambda x, y: x | y, (1 << j for j in g))
+
+            # the bitmap set of all tensors is computed as it is needed to
+            # compute set differences: s1 - s2 transforms into
+            # s1 & (all_tensors ^ s2)
+            all_tensors = (1 << len(inputs)) - 1
+
+            # try to find contraction with cost <= cost_cap and increase
+            # cost_cap successively if no such contraction is found;
+            # this is a major performance improvement; start with product of
+            # output index dimensions as initial cost_cap
+            cost_cap = helpers.compute_size_by_dict(
+                set.union(*(inputs[j] for j in _bitmapset_indices(g))) & output,
+                size_dict
+            )
+
+            while len(x[-1]) == 0:
+                for n in range(2, len(x[1]) + 1):
+                    xn = x[n]
+
+                    # try to combine solutions from x[m] and x[n-m]
+                    for m in range(1, n // 2 + 1):
+                        for s1, (i1, cost1, cntrct1) in x[m].items():
+                            for s2, (i2, cost2, cntrct2) in x[n-m].items():
+
+                                # only if s1 and s2 are disjoint
+                                if s1 & s2 == 0:
+
+                                    # avoid e.g. s1={0}, s2={1} and s1={1}, s2={0}
+                                    if m != n - m or s1 < s2:
+
+                                        i1_cut_i2_wo_output = (i1 & i2) - output
+
+                                        # ignore outer products:
+                                        if len(i1_cut_i2_wo_output) > 0:
+
+                                            i1_union_i2 = i1 | i2
+                                            cost = cost1 + cost2 + helpers.compute_size_by_dict(i1_union_i2, size_dict)
+                                            if cost <= cost_cap:
+                                                s = s1 | s2
+                                                if s not in xn or cost < xn[s][1]:
+                                                    # set of remaining tensors (=g-s)
+                                                    r = g & (all_tensors ^ s)
+
+                                                    # indices of remaining indices:
+                                                    i_r = (set.union(*(inputs[j] for j in _bitmapset_indices(r)))
+                                                           if r != 0 else set())
+
+                                                    # contraction indices:
+                                                    i_cntrct = i1_cut_i2_wo_output - i_r
+
+                                                    i = i1_union_i2 - i_cntrct
+                                                    mem = helpers.compute_size_by_dict(i, size_dict)
+                                                    if memory_limit is None or mem <= memory_limit:
+                                                        xn[s] = (i, cost, (cntrct1, cntrct2))
+
+                # increase cost cap for next iteration:
+                cost_cap = min(size_dict) * cost_cap
+
+            i, cost, contraction = list(x[-1].values())[0]
+            subgraph_contractions.append(contraction)
+            subgraph_contractions_size.append(helpers.compute_size_by_dict(i, size_dict))
+
+        # sort the subgraph contractions by the size of the subgraphs in
+        # ascending order (will give the cheapest contractions); note that
+        # outer products should be performed pairwise (to use BLAS functions)
+        subgraph_contractions = [
+            subgraph_contractions[j]
+            for j in np.argsort(subgraph_contractions_size)
+        ]
+
+        # build the final contraction tree
+        tree = functools.reduce(lambda x, y: (x, y), subgraph_contractions)
+
+        return _tree_to_sequence(tree)
+
+
+def dynamic_programming(inputs, output, size_dict, memory_limit=None):
+    optimizer = DynamicProgrammingOptimizer()
+    return optimizer(inputs, output, size_dict, memory_limit)
+
+
 _AUTO_CHOICES = {}
 for i in range(1, 5):
     _AUTO_CHOICES[i] = optimal
@@ -680,6 +1001,8 @@ _PATH_OPTIONS = {
     'greedy': greedy,
     'eager': greedy,
     'opportunistic': greedy,
+    'dp': dynamic_programming,
+    'dynamic-programming': dynamic_programming
 }
 
 


Test patch:
diff --git a/opt_einsum/tests/test_contract.py b/opt_einsum/tests/test_contract.py
index 643f724..bc0d4c6 100644
--- a/opt_einsum/tests/test_contract.py
+++ b/opt_einsum/tests/test_contract.py
@@ -97,7 +97,7 @@ tests = [
 ]
 
 
-all_optimizers = ['optimal', 'branch-all', 'branch-2', 'branch-1', 'greedy', 'random-greedy']
+all_optimizers = ['optimal', 'branch-all', 'branch-2', 'branch-1', 'greedy', 'random-greedy', 'dp']
 
 
 @pytest.mark.parametrize("string", tests)
diff --git a/opt_einsum/tests/test_paths.py b/opt_einsum/tests/test_paths.py
index ae0cc1f..dfddf27 100644
--- a/opt_einsum/tests/test_paths.py
+++ b/opt_einsum/tests/test_paths.py
@@ -26,23 +26,29 @@ explicit_path_tests = {
     }),
 }
 
+# note that these tests have no unique solution due to the chosen dimensions
 path_edge_tests = [
     ['greedy', 'eb,cb,fb->cef', ((0, 2), (0, 1))],
     ['branch-all', 'eb,cb,fb->cef', ((0, 2), (0, 1))],
     ['branch-2', 'eb,cb,fb->cef', ((0, 2), (0, 1))],
     ['optimal', 'eb,cb,fb->cef', ((0, 2), (0, 1))],
+    ['dp', 'eb,cb,fb->cef', ((1, 2), (0, 1))],
     ['greedy', 'dd,fb,be,cdb->cef', ((0, 3), (0, 1), (0, 1))],
     ['branch-all', 'dd,fb,be,cdb->cef', ((0, 3), (0, 1), (0, 1))],
     ['branch-2', 'dd,fb,be,cdb->cef', ((0, 3), (0, 1), (0, 1))],
     ['optimal', 'dd,fb,be,cdb->cef', ((0, 3), (0, 1), (0, 1))],
+    ['optimal', 'dd,fb,be,cdb->cef', ((0, 3), (0, 1), (0, 1))],
+    ['dp', 'dd,fb,be,cdb->cef', ((0, 3), (0, 2), (0, 1))],
     ['greedy', 'bca,cdb,dbf,afc->', ((1, 2), (0, 2), (0, 1))],
     ['branch-all', 'bca,cdb,dbf,afc->', ((1, 2), (0, 2), (0, 1))],
     ['branch-2', 'bca,cdb,dbf,afc->', ((1, 2), (0, 2), (0, 1))],
     ['optimal', 'bca,cdb,dbf,afc->', ((1, 2), (0, 2), (0, 1))],
+    ['dp', 'bca,cdb,dbf,afc->', ((1, 2), (1, 2), (0, 1))],
     ['greedy', 'dcc,fce,ea,dbf->ab', ((1, 2), (0, 1), (0, 1))],
     ['branch-all', 'dcc,fce,ea,dbf->ab', ((1, 2), (0, 2), (0, 1))],
     ['branch-2', 'dcc,fce,ea,dbf->ab', ((1, 2), (0, 2), (0, 1))],
     ['optimal', 'dcc,fce,ea,dbf->ab', ((1, 2), (0, 2), (0, 1))],
+    ['dp', 'dcc,fce,ea,dbf->ab', ((1, 2), (0, 2), (0, 1))],
 ]
 
 
@@ -191,7 +197,7 @@ def test_greedy_edge_cases():
     assert check_path(path, [(0, 1), (0, 2), (0, 1)])
 
 
-@pytest.mark.parametrize("optimize", ['greedy', 'branch-2', 'branch-all', 'optimal'])
+@pytest.mark.parametrize("optimize", ['greedy', 'branch-2', 'branch-all', 'optimal', 'dp'])
 def test_can_optimize_outer_products(optimize):
     a, b, c = [np.random.randn(10, 10) for _ in range(3)]
     d = np.random.randn(10, 2)