Ellipsis Broadcasting at numpy.einsum

Question

I am having trouble understanding why the following does not work:

I have an array prefactor that can be three-dimensional or six-dimensional. I have dipole matrices that have four dimensions. The first three dimensions of the dipoles correspond to the last three dimensions of the pre-factor.

Since I do not know the form of the prefactor, I use Ellipsis to account for the three optional dimensions in the prefactor:

numpy.einsum('...lmn,lmno->...o', prefactor, dipoles) 

(In the example here, prefactor.shape (1, 1, 1, 160, 160, 128) and dipoles.shape (160, 160, 128, 3). When I execute, I get the error:

operand 1 was not large enough to fit the broadcast and could not be expanded, since the indices of the Einstein indices were indicated both at the beginning and at the end

This works, however, when I add an ellipsis to the second term:

 numpy.einsum('...lmn,...lmno->...o', prefactor, dipoles) 

I just don’t understand why, because there should not be a need for an ellipse. Does anyone know what is going on here?

The same question was asked at http://comments.gmane.org/gmane.comp.python.numeric.general/53705 , but there is no satisfactory answer yet.

Answer 1

There is a github problem for this problem:

https://github.com/numpy/numpy/issues/2455 improving index notation in einsum (Tra # 1862)

Error Case:

 einsum('ij...,j->ij...',A,B) 

The current job requires an (empty) ellipse:

einsum ('IJ ..., J ...-> IJ ...', A, B)

It seems that einsum repeated through the string and ops arguments several times, identifying the indices and broadcast types (right, left, middle, none) and op. In doing so, he builds a numpy.nditer . It is when building op_axes for nditer that einsum is causing this error. I don’t know if the test criteria are too complicated ( ibroadcast >= ndim ), or if one more step is necessary for this argument to build the correct op_axes .

https://github.com/numpy/numpy/issues/2619 shows how nditer can be used to replicate einsum behavior. Working with this, I can reproduce your calculations this way:

 prefactor = np.random.random((1, 1, 1, 160, 160, 128)) dipoles = np.random.random((160, 160, 128, 3)) x = numpy.einsum('...lmn,...lmno->...o', prefactor, dipoles) #numpy.einsum('...lmn,lmno->...o', prefactor, dipoles) # not work op_axes = [[0,1,2,3,4,5,-1], [-1,-1,-1,0,1,2,3], [0,1,2,-1,-1,-1,3]] flags = ['reduce_ok','buffered', 'external_loop', 'delay_bufalloc', 'grow_inner'] op_flags = [['readonly']]*nops + [['allocate','readwrite']] it = np.nditer([prefactor,dipoles,None], flags, op_flags, op_axes=op_axes) it.operands[nops][...] = 0 it.reset() #it.debug_print() for (x,y,w) in it: w[...] += x*y print "\nnditer usage:" print it.operands[nops] # == x print it.operands[nops].shape # (1, 1, 1, 3) 

The op_axes line indicates that einsum outputting from '...lmn,...lmno->...o' .

---

I am studying this problem at https://github.com/hpaulj/numpy-einsum .

There I have einsum_py.py that emulates np.einsum with Python code. The part related to this problem is parse_subscripts() , and in particular prepare_op_axes() . It seems that for the correct creation of op_axes only the iteration BROADCAST_RIGHT is required (starting from the end), regardless of where the ellipses are in the indices. It also removes the error message underlying this problem.

The einsum.c.src file in this repository has this change and compiles correctly with the current main distribution (just replace the file and assembly). It is well versed in relation to test_einsum.py , as well as examples from this problem.

I sent a transfer request for this change.