too many subscripts in einsum

[punkdit]

Here is the bottom of my stack trace:

  File "/home/simon/site-packages/opt_einsum/contract.py", line 731, in __call__
    return self._contract(ops, out, backend, evaluate_constants=evaluate_constants)
  File "/home/simon/site-packages/opt_einsum/contract.py", line 667, in _contract
    **self.einsum_kwargs)
  File "/home/simon/site-packages/opt_einsum/contract.py", line 563, in _core_contract
    new_view = _einsum(einsum_str, *tmp_operands, backend=backend, **einsum_kwargs)
  File "/home/simon/site-packages/opt_einsum/sharing.py", line 161, in cached_einsum
    return einsum(*args, **kwargs)
  File "/home/simon/site-packages/opt_einsum/contract.py", line 360, in _einsum
    return fn(einsum_str, *operands, **kwargs)
ValueError: Internal error while evaluating `ContractExpression`. Note that few checks are performed - the number and rank of the array arguments must match the original expression. The internal error was: '('too many subscripts in einsum',)'

It looks like numpy-einsum is barfing here:

def iterlabels(labelstr, output_labels, ndim_output):
    """
     * Set up the labels for the iterator (output + combined labels).
     * Can just share the output_labels memory, because iter_labels
     * is output_labels with some more labels appended.
    """
    iter_labels = output_labels[:]
    ndim_iter = ndim_output
    for label in range(labelstr.min_label, labelstr.max_label+1):
        if labelstr.counts.get(label, 0) > 0 and label not in output_labels:
            if ndim_iter > 128:
                raise ValueError('too many subscripts in einsum')
            iter_labels.append(label)
            ndim_iter += 1
    # may have 0 in iter_labels
    return iter_labels

This is strange because I get the path_info from contract_path and it tells me
that largest_intermediate is of size 134217728, which seems doable.

I'm wondering how einsum ended up with more that 128 labels (i'm assuming these are tensor indices?) Perhaps there are a bunch of trivial indices being summed over...

Any suggestions ?

[dgasmith]

An intermediate size of ~1GB (doubles) is certainly fine :), but the issue here is the number of indices. Assuming that each index is only of size two leads us to an impractical computational cost.

This is a bug, can you post a reproducible example?

[jcmgray]

Any chance of posting the input contraction, even in a gist or something? If the size is 2**27 with 128 indices, many must be size 1 in which case some squeezing might be in order.

As another check, opt_einsum prefers einsum for some operations even though it is not necessary, you might try setting

from opt_einsum.backends.dispatch import _has_einsum
_has_einsum['numpy'] = False

to force it to use tensordot (though if the contraction requires einsum, for e.g. 'ii->i', then there is still a problem).

[punkdit]

Here it is:

import numpy
import opt_einsum as oe

shapes = [(2,), (2,), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2,), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2,), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2), (2, 2),
    (2, 2), (2, 2), (2, 2), (2, 2)]

linkss = [[10], [0], [0, 11], [0, 1], [1, 12], [1, 2], [2, 13],
    [2, 3], [3, 14], [3, 4], [4, 15], [4, 5], [5, 16], [5,
    6], [6, 17], [6, 7], [7, 18], [7, 8], [8, 19], [8, 9],
    [9, 20], [9], [10, 21], [10, 11], [11, 22], [11, 12],
    [12, 23], [12, 13], [13, 24], [13, 14], [14, 25], [14,
    15], [15, 26], [15, 16], [16, 27], [16, 17], [17, 28],
    [17, 18], [18, 29], [18, 19], [19, 30], [19, 20], [20,
    31], [10, 20], [21, 32], [21, 22], [22, 33], [22, 23],
    [23, 34], [23, 24], [24, 35], [24, 25], [25, 36], [25,
    26], [26, 37], [26, 27], [27, 38], [27, 28], [28, 39],
    [28, 29], [29, 40], [29, 30], [30, 41], [30, 31], [31,
    42], [21, 31], [32, 43], [32, 33], [33, 44], [33, 34],
    [34, 45], [34, 35], [35, 46], [35, 36], [36, 47], [36,
    37], [37, 48], [37, 38], [38, 49], [38, 39], [39, 50],
    [39, 40], [40, 51], [40, 41], [41, 52], [41, 42], [42,
    53], [32, 42], [43, 54], [43, 44], [44, 55], [44, 45],
    [45, 56], [45, 46], [46, 57], [46, 47], [47, 58], [47,
    48], [48, 59], [48, 49], [49, 60], [49, 50], [50, 61],
    [50, 51], [51, 62], [51, 52], [52, 63], [52, 53], [53,
    64], [43, 53], [54, 65], [54, 55], [55, 66], [55, 56],
    [56, 67], [56, 57], [57, 68], [57, 58], [58, 69], [58,
    59], [59, 70], [59, 60], [60, 71], [60, 61], [61, 72],
    [61, 62], [62, 73], [62, 63], [63, 74], [63, 64], [64,
    75], [54, 64], [65, 76], [65, 66], [66, 77], [66, 67],
    [67, 78], [67, 68], [68, 79], [68, 69], [69, 80], [69,
    70], [70, 81], [70, 71], [71, 82], [71, 72], [72, 83],
    [72, 73], [73, 84], [73, 74], [74, 85], [74, 75], [75,
    86], [65, 75], [76, 87], [76, 77], [77, 88], [77, 78],
    [78, 89], [78, 79], [79, 90], [79, 80], [80, 91], [80,
    81], [81, 92], [81, 82], [82, 93], [82, 83], [83, 94],
    [83, 84], [84, 95], [84, 85], [85, 96], [85, 86], [86,
    97], [76, 86], [87, 98], [87, 88], [88, 99], [88, 89],
    [89, 100], [89, 90], [90, 101], [90, 91], [91, 102],
    [91, 92], [92, 103], [92, 93], [93, 104], [93, 94], [94,
    105], [94, 95], [95, 106], [95, 96], [96, 107], [96,
    97], [97, 108], [87, 97], [98, 109], [98, 99], [99, 110],
    [99, 100], [100, 111], [100, 101], [101, 112], [101,
    102], [102, 113], [102, 103], [103, 114], [103, 104],
    [104, 115], [104, 105], [105, 116], [105, 106], [106,
    117], [106, 107], [107, 118], [107, 108], [108, 119],
    [98, 108], [109], [109, 110], [0, 110], [110, 111], [1,
    111], [111, 112], [2, 112], [112, 113], [3, 113], [113,
    114], [4, 114], [114, 115], [5, 115], [115, 116], [6,
    116], [116, 117], [7, 117], [117, 118], [8, 118], [118,
    119], [9, 119], [109, 119]]


As = [numpy.zeros(shape) for shape in shapes]
n = len(As)

args = []
for i in range(n):
    args.append(As[i])
    args.append(linkss[i])

v = oe.contract(*args)

[jcmgray]

So the issue seems to be that numpy.einsum doesn't seem to be able to handle more than 26 indices on any input. This is a proxy for the step that breaks:

import numpy as np

A = np.random.rand(*[2]*27)
B = np.random.rand(*[2]*14)
eq = 'jbcdreoGEklmxyvzABCqnpFDwtu,hibqpgmnlafkos->jbcdreGExyvzABCFDwtuhigafs'

np.einsum(eq, A, B)

Unfortunately, this step has indices that appear more than twice ('b') so cannot be forced to use tensordot.

I think your options are:

1. Use the 'random-greedy' optimizer, this produces a much better contraction (2**22 max size).
2. Raise an issue with numpy einsum proper.

[dgasmith]

Odd, I had thought it would be 52 indices (a-zA-Z). I need to poke through their C code for something else, I might try to see what I can do.

[bbbbbbbbba]

It isn't the number of indices in a single input, it's the total number of indices involved (more precisely, the number of indices that numpy needs to iterate over). The relevant C code in numpy is:

https://github.com/numpy/numpy/blob/4cba2d91e1546872d29af6b25ad35947f27e03ac/numpy/core/src/multiarray/einsum.c.src#L947-L964

    /*
     * Set up the labels for the iterator (output + combined labels).
     * Can just share the output_labels memory, because iter_labels
     * is output_labels with some more labels appended.
     */
    iter_labels = output_labels;
    ndim_iter = ndim_output;
    for (label = min_label; label <= max_label; ++label) {
        if (label_counts[label] > 0 &&
                memchr(output_labels, label, ndim_output) == NULL) {
            if (ndim_iter >= NPY_MAXDIMS) {
                PyErr_SetString(PyExc_ValueError,
                            "too many subscripts in einsum");
                goto fail;
            }
            iter_labels[ndim_iter++] = label;
        }
    }

https://github.com/numpy/numpy/blob/8281acdb2900d693cf2f8777c3707e502bae89d9/numpy/core/include/numpy/ndarraytypes.h#L40

#define NPY_MAXDIMS 32

[dgasmith]

Thanks for dropping a comment, this makes perfect sense. A dim size of anything but basically 2 at this size would not be possible to compute.

[dgasmith]

@jcmgray reported as a NumPy issue, not a lot we can do here.