16

Consider the following exercise in Numpy array broadcasting.

import numpy as np
v = np.array([[1.0, 2.0]]).T # column array

A2 = np.random.randn(2,10) # 2D array
A3 = np.random.randn(2,10,10) # 3D

v * A2 # works great

# causes error: 
v * A3 # error

I know the Numpy rules for broadcasting, and I'm familiar with bsxfun functionality in Matlab. I understand why attempting to broadcast a (2,1) array into a (2,N,N) array fails, and that I have to reshape the (2,1) array into a (2,1,1) array before this broadcasting goes through.

My question is: is there any way to tell Python to automatically pad the dimensionality of an array when it attempts to broadcast, without me having to specifically tell it the necessary dimension?

I don't want to explicitly couple the (2,1) vector with the multidimensional array it's going to be broadcast against---otherwise I could do something stupid and absurdly ugly like mult_v_A = lambda v,A: v.reshape([v.size] + [1]*(A.ndim-1)) * A. I don't know ahead of time if the "A" array will be 2D or 3D or N-D.

Matlab's bsxfun broadcasting functionality implicitly pads the dimensions as needed, so I'm hoping there's something I could do in Python.

Improve this question
6
  • 4
    My view is that this is a feature and not a bug. For example, suppose you had a 2-by-1 column vector v and then you had a 2-by-2-by-10 ndarray. Do you want to reshape v to have shape (2,1,1) or shape (1,2,1)? If you just pad the dimensions, it could be ambiguous to the user. Forcing an explicit reshape is a better general procedure, and leave it to the user to write a special function to perform the reshaping automatically if the user has a fixed convention. But it's not good to make a global numpy dimension-padder that forces a convention upon you. It would be too easy to misuse. Commented Jul 10, 2013 at 17:48
  • -1 @EMS. This ambiguity is readily solved by specifying that the first non-singleton dimension will be used in the broadcast. This attitude of "This way it's better" is totally inappropriate for systems used by professional programmers and applied mathematicians---this is a flaw in terms of illegibility and complectedness, not a feature. Commented Jul 10, 2013 at 18:04
  • Also, though the questions start out differently, this is a possible duplicate of: < stackoverflow.com/questions/15296944/… >. I think numpy.apply_along_axis and numpy.tensordot in particular seem more than sufficient for this kind of thing, while still leaving the nice property that the programmer must make some explicit reference to the way that ambiguous dimension changes should occur for broadcasting. Commented Jul 10, 2013 at 18:40
  • I'd expect apply_along_axis to totally destroy performance by de-vectorizing the array operation. Also, broadcasting is much more useful than just multiplication via tensordot. Commented Jul 10, 2013 at 18:47
  • @EMS, I expected you to be a professional user---it's hard to see foibles in one's daily tools as foibles; instead we like to praise them. As unutbu's answer illustrates, Numpy has already chosen to resolve such ambiguities, using simply the reverse of my suggestion: Numpy is perfectly happy broadcasting (10, 2, 2) and (2, 1) arrays, and broadcasts over the last non-singleton dimension because it works backwards along the dimension tuple. Commented Jul 10, 2013 at 18:56

2 Answers 2

Reset to default
11

It's ugly, but this will work:

(v.T * A3.T).T

If you don't give it any arguments, transposing reverses the shape tuple, so you can now rely on the broadcasting rules to do their magic. The last transpose returns everything to the right order.

Improve this answer
Sign up to request clarification or add additional context in comments.

1 Comment

It is actually less ugly than what I had after studying @unutbu's answer, which was np.swapaxes(v.T * np.swapaxes(A3,0,-1), 0, -1)! Knowing that transpose reverses all dimensions makes this much neater. At least users will see these transposes and hopefully think "Ok there's broadcasting going on here, with just a tiny touch of funkiness to make the dimensions line up, ok, I'm ok". Thanks chief.
9

NumPy broadcasting adds additional axes on the left.

So if you arrange your arrays so the shared axes are on the right and the broadcastable axes are on the left, then you can use broadcasting with no problem:

import numpy as np
v = np.array([[1.0, 2.0]])  # shape (1, 2)

A2 = np.random.randn(10,2) # shape (10, 2)
A3 = np.random.randn(10,10,2) # shape (10, 10, 2)

v * A2  # shape (10, 2)

v * A3 # shape (10, 10, 2)
Improve this answer

4 Comments

This at least is a half-way solution. The established convention in my application is unfortunately that vectors like v are column arrays. Assuming I didn't want to reorganize the code (and my users) to expect row arrays, I am debating whether a function that swaps dimensions to do a broadcast multiply is really any less complected than a reshape operation that intertwines A and v in my original post.
Is there some reason why Numpy adds dimensions to the left, but not on the right? It seems arbitrary but perhaps there's some important rationale related to memory or algebra?
Possibly answering my own (sub-)question, is it because broadcasting starts at arrays' trailing dimensions and works leftward, and keeps adding dimensions as needed, as long as the rules are met? Interesting!
NumPy arrays are in C-contiguous order (last-index varies the fastest) by default. So adding additional axes on the left is essentially costless -- only the strides and shape attributes have to change, while the underlying data can remain the same. If NumPy arrays were in F-contiguous order by default, then broadcasting would work most naturally by adding axes on the right.

Your Answer

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.