1

The code below is meant to conduct a linear coordinate transformation on a set of 3d coordinates. The transformation matrix is A, and the array containing the coordinates is x. The zeroth axis of x runs over the dimensions x, y, z. It can have any arbitrary shape beyond that.

Here's my attempt:

A = np.random.random((3, 3))
x = np.random.random((3, 4, 2))

x_prime = np.einsum('ij,j...->i...', A, x)

The output is:

    x_prime = np.einsum('ij,j...->i...', A, x)

ValueError: operand 0 did not have enough dimensions
to match the broadcasting, and couldn't be extended
because einstein sum subscripts were specified at both
the start and end

If I specify the additional subscripts in x explicitly, the error goes away. In other words, the following works:

x_prime = np.einsum('ij,jkl->ikl', A, x)

I'd like x to be able to have any arbitrary number of axes after the zeroth axis, so the workaround I give about is not optimal. I'm actually not sure why the first einsum example is not working. I'm using numpy 1.6.1. Is this a bug, or am I misunderstanding the documentation?

apdnu
  • 6,717
  • 5
  • 24
  • 24
  • FWIW your code seems to work as-is in numpy 1.9.0.dev-b785070. – DSM Feb 13 '14 at 18:31
  • I've noticed at least one other bug with einsum in my version of numpy (incorrect results with three operands), so I wouldn't be surprised if this is a second bug. I think einsum was introduced in v1.6.0. – apdnu Feb 13 '14 at 18:39
  • `x_prime = np.einsum('...ij,j...->i...', A, x)` should also work. This issue was raised in `http://stackoverflow.com/questions/16591696/ellipsis-broadcasting-in-numpy-einsum/`, and corrected in the latest code. – hpaulj Feb 13 '14 at 19:05

1 Answers1

3

Yep, it's a bug. It was fixed in this pull request: https://github.com/numpy/numpy/pull/4099

This was only merged a month ago, so it'll be a while before it makes it to a stable release.

EDIT: As @hpaulj mentions in the comment, you can work around this limitation by adding an ellipsis even when all indices are specified:

np.einsum('...ij,j...->i...', A, x)
perimosocordiae
  • 17,287
  • 14
  • 60
  • 76
  • If you add the suggestion from @hpaulj's comment, this answer would be spot on. Other people will probably have the same problem for some time to come, since it will take a while for the next stable version of numpy to get included in the major distros (Ubuntu 14.04 LTS probably won't include it, CentOS takes ages to release new versions, etc.). – apdnu Feb 13 '14 at 22:50