title: Multidimensional Broadcasting in Python / NumPy - or inverse of `numpy.squeeze()`
body: <p>What would be the best way of broadcasting two arrays together when a simple call to <code>np.broadcast_to()</code> would fail?</p>
  
  <p>Consider the following example:</p>
  
  <pre><code>import numpy as np
  
  arr1 = np.arange(2 * 3 * 4 * 5 * 6).reshape((2, 3, 4, 5, 6))
  arr2 = np.arange(3 * 5).reshape((3, 5))
  
  arr1 + arr2
  # ValueError: operands could not be broadcast together with shapes (2,3,4,5,6) (3,5)
  
  arr2_ = np.broadcast_to(arr2, arr1.shape)
  # ValueError: operands could not be broadcast together with remapped shapes
  
  arr2_ = arr2.reshape((1, 3, 1, 5, 1))
  arr1 + arr2
  # now this works because the singletons trigger the automatic broadcast
  </code></pre>
  
  <p>This only work if I manually select a shape for which automatic broadcasting is going to work.
  What would be the most efficient way of doing this automatically?
  Is there an alternative way other than reshape on a cleverly constructed <em>broadcastable</em> shape?</p>
  
  <p>Note the relation to <code>np.squeeze()</code>: this would perform the inverse operation by removing singletons. So what I need is some sort of <code>np.squeeze()</code> inverse.
  The official <a href="https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.squeeze.html" rel="nofollow noreferrer">documentation</a> (as of <em>NumPy 1.13.0</em> suggests that the inverse of <code>np.squeeze()</code> is <code>np.expand_dim()</code>, but this is not nearly as flexible as I'd need it to be, and actually <code>np.expand_dim()</code> is roughly equivalent to <code>np.reshape(array, shape + (1,))</code> or <code>array[:, None]</code>.</p>
  
  <p>This issue is also related to the <code>keepdims</code> keyword accepted by e.g. <code>sum</code>:</p>
  
  <pre><code>import numpy as np
  
  arr1 = np.arange(2 * 3 * 4 * 5 * 6).reshape((2, 3, 4, 5, 6))
  
  # not using `keepdims`
  arr2 = np.sum(arr1, (0, 2, 4))
  arr2.shape
  # : (3, 5)
  
  arr1 + arr2
  # ValueError: operands could not be broadcast together with shapes (2,3,4,5,6) (3,5)
  
  
  # now using `keepdims`
  arr2 = np.sum(arr1, (0, 2, 4), keepdims=True)
  arr2.shape
  # : (1, 3, 1, 5, 1)
  
  arr1 + arr2
  # now this works because it has the correct shape
  </code></pre>
  
  <p><strong>EDIT</strong>: Obviously, in cases where <code>np.newaxis</code> or <code>keepdims</code> mechanisms are an appropriate choice, there would be no need for a <code>unsqueeze()</code> function.</p>
  
  <p>Yet, there are use-cases where none of these is an option.</p>
  
  <p>For example, consider the case of the weighted average as implemented in <code>numpy.average()</code> over an arbitrary number of dimensions specified by <code>axis</code>.
  Right now the <code>weights</code> parameter must have the same shape as the input.
  However, <code>weights</code> there is no need specify the weights over the non-reduced dimensions as they are just repeating and the NumPy's broadcasting mechanism would appropriately take care of them.</p>
  
  <p>So if we would like to have such a functionality, we would need to code something like (where some consistency checks are just omitted for simplicity):</p>
  
  <pre><code>def weighted_average(arr, weights=None, axis=None):
      if weights is not None and weights.shape != arr.shape:
          weights = unsqueeze(weights, ...)
          weights = np.zeros_like(arr) + weights
      result = np.sum(arr * weights, axis=axis)
      result /= np.sum(weights, axis=axis)
      return result
  </code></pre>
  
  <p>or, equivalently:</p>
  
  <pre><code>def weighted_average(arr, weights=None, axis=None):
      if weights is not None and weights.shape != arr.shape:
          weights = unsqueeze(weights, ...)
          weights = np.zeros_like(arr) + weights
      return np.average(arr, weights, axis)
  </code></pre>
  
  <p>In either of the two, it is not possible to replace <code>unsqueeze()</code> with <code>weights[:, np.newaxis]</code>-like statements because we do not know beforehand where the new axis will be needed, nor we can use the <code>keepdims</code> feature of <code>sum</code> because the code will fail at <code>arr * weights</code>.</p>
  
  <p>This case could be relatively nicely handled if <a href="https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.expand_dims.html" rel="nofollow noreferrer"><code>np.expand_dims()</code></a> would support an iterable of ints for its <code>axis</code> parameter, but as of NumPy 1.13.0 does not.</p>