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?

3
  • FWIW your code seems to work as-is in numpy 1.9.0.dev-b785070.
    – DSM
    Feb 13, 2014 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, 2014 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, 2014 at 19:05

1 Answer 1

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)
1
  • 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, 2014 at 22:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.