1

What would be the best way of broadcasting two arrays together when a simple call to np.broadcast_to() would fail?

Consider the following example:

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

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 broadcastable shape?

Note the relation to np.squeeze(): this would perform the inverse operation by removing singletons. So what I need is some sort of np.squeeze() inverse. The official documentation (as of NumPy 1.13.0 suggests that the inverse of np.squeeze() is np.expand_dim(), but this is not nearly as flexible as I'd need it to be, and actually np.expand_dim() is roughly equivalent to np.reshape(array, shape + (1,)) or array[:, None].

This issue is also related to the keepdims keyword accepted by e.g. sum:

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

EDIT: Obviously, in cases where np.newaxis or keepdims mechanisms are an appropriate choice, there would be no need for a unsqueeze() function.

Yet, there are use-cases where none of these is an option.

For example, consider the case of the weighted average as implemented in numpy.average() over an arbitrary number of dimensions specified by axis. Right now the weights parameter must have the same shape as the input. However, weights 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.

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):

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

or, equivalently:

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)

In either of the two, it is not possible to replace unsqueeze() with weights[:, np.newaxis]-like statements because we do not know beforehand where the new axis will be needed, nor we can use the keepdims feature of sum because the code will fail at arr * weights.

This case could be relatively nicely handled if np.expand_dims() would support an iterable of ints for its axis parameter, but as of NumPy 1.13.0 does not.

norok2
  • 25,683
  • 4
  • 73
  • 99
  • 2
    This is horribly unsafe (and `np.squeeze` is pretty unsafe too, but this is worse). If any inputs happen to have multiple dimensions of the same length, you have no way to tell how to align the shapes of the two arrays. It's much better to insert the new axes in the correct positions with `np.newaxis`, or use mechanisms like `keepdims` to ensure the shapes are aligned in the first place. – user2357112 Sep 26 '17 at 17:45
  • Your concerns are well taken, but I clearly see some example where neither is appropriate. I have added an example where neither `np.squeeze` nor `keepdims` could be used. – norok2 Sep 26 '17 at 18:26
  • That's still a case for `np.newaxis`. You might have to construct the indexing tuple manually, or you could use `reshape` and construct the shape tuple manually, but it's definitely not a job for `unsqueeze`. – user2357112 Sep 26 '17 at 18:30
  • (Or, of course, you could have the caller align the shapes. That's probably better, since it's not actually clear that the weights should be repeated over the non-reduced dimensions instead of following the normal broadcasting rules.) – user2357112 Sep 26 '17 at 18:31
  • what is the difference between constructing the tuple "manually" and having a function doing that for you? Or, alternatively, what would be the rationale for the manual construction? I am really struggling to see that. – norok2 Sep 26 '17 at 18:32
  • 2
    Constructing it manually, you can determine where the new axes have to go from information other than the existing shapes. For example, `weighted_average` could use the `axis` argument, or the caller could use the information they have to align the shapes. This is much safer than guessing based on the lengths of the dimensions. – user2357112 Sep 26 '17 at 18:34
  • But that would mean that some `axis` argument fine-tuning the `unsqueeze` behavior is desireable. This does not mean that the approach of dynamically adding singletons is *per se* a bad idea. – norok2 Sep 26 '17 at 18:41
  • @user2357122 What is unsafe about using `np.squeeze`? – user2699 Sep 26 '17 at 18:41
  • 2
    @user2699: If more dimensions than you expected have length 1, you'll end up with the wrong output shape. – user2357112 Sep 26 '17 at 18:43
  • @user2357112, good point, but I don't see how that's more dangerous than any other operation would be when the dimensions of the input aren't what you thought. – user2699 Sep 26 '17 at 20:20
  • @user2699: 1 is an entirely normal, valid length for a dimension to have, and most code will handle length-1 dimensions just fine even if you don't explicitly think about it. You don't have to special-case a data set with only one data point, or observations with only one measured variable, even if your code was written to handle large inputs with many measured variables... unless you used `np.squeeze`. `np.squeeze` will squeeze any length-1 dimensions without any idea why they were length-1 or whether you wanted to keep some of them. – user2357112 Sep 26 '17 at 20:30
  • `expand_dim` takes an `axis` parameter, so does what you do, but only one singleton at a time. – hpaulj Sep 27 '17 at 08:10
  • That's the whole point, the `axis` parameter of `expand_dim()` only supports `int`. I have uploaded my answer so that now `unsqueeze()` handles an `axis` parameters which can be an iterable of ints. – norok2 Sep 27 '17 at 12:16

1 Answers1

3

My way of achieving this is by defining the following unsqueezing() function to handle cases where this can be done automatically and giving a warning when the inputs could be ambiguous (e.g. when some source elements of the source shape may match multiple elements of the target shape):

def unsqueezing(
        source_shape,
        target_shape):
    """
    Generate a broadcasting-compatible shape.

    The resulting shape contains *singletons* (i.e. `1`) for non-matching dims.
    Assumes all elements of the source shape are contained in the target shape
    (excepts for singletons) in the correct order.

    Warning! The generated shape may not be unique if some of the elements
    from the source shape are present multiple timesin the target shape.

    Args:
        source_shape (Sequence): The source shape.
        target_shape (Sequence): The target shape.

    Returns:
        shape (tuple): The broadcast-safe shape.

    Raises:
        ValueError: if elements of `source_shape` are not in `target_shape`.

    Examples:
        For non-repeating elements, `unsqueezing()` is always well-defined:

        >>> unsqueezing((2, 3), (2, 3, 4))
        (2, 3, 1)
        >>> unsqueezing((3, 4), (2, 3, 4))
        (1, 3, 4)
        >>> unsqueezing((3, 5), (2, 3, 4, 5, 6))
        (1, 3, 1, 5, 1)
        >>> unsqueezing((1, 3, 5, 1), (2, 3, 4, 5, 6))
        (1, 3, 1, 5, 1)

        If there is nothing to unsqueeze, the `source_shape` is returned:

        >>> unsqueezing((1, 3, 1, 5, 1), (2, 3, 4, 5, 6))
        (1, 3, 1, 5, 1)
        >>> unsqueezing((2, 3), (2, 3))
        (2, 3)

        If some elements in `source_shape` are repeating in `target_shape`,
        a user warning will be issued:

        >>> unsqueezing((2, 2), (2, 2, 2, 2, 2))
        (2, 2, 1, 1, 1)
        >>> unsqueezing((2, 2), (2, 3, 2, 2, 2))
        (2, 1, 2, 1, 1)

        If some elements of `source_shape` are not presente in `target_shape`,
        an error is raised.

        >>> unsqueezing((2, 3), (2, 2, 2, 2, 2))
        Traceback (most recent call last):
          ...
        ValueError: Target shape must contain all source shape elements\
 (in correct order). (2, 3) -> (2, 2, 2, 2, 2)
        >>> unsqueezing((5, 3), (2, 3, 4, 5, 6))
        Traceback (most recent call last):
          ...
        ValueError: Target shape must contain all source shape elements\
 (in correct order). (5, 3) -> (2, 3, 4, 5, 6)

    """
    shape = []
    j = 0
    for i, dim in enumerate(target_shape):
        if j < len(source_shape):
            shape.append(dim if dim == source_shape[j] else 1)
            if i + 1 < len(target_shape) and dim == source_shape[j] \
                    and dim != 1 and dim in target_shape[i + 1:]:
                text = ('Multiple positions (e.g. {} and {})'
                        ' for source shape element {}.'.format(
                    i, target_shape[i + 1:].index(dim) + (i + 1), dim))
                warnings.warn(text)
            if dim == source_shape[j] or source_shape[j] == 1:
                j += 1
        else:
            shape.append(1)
    if j < len(source_shape):
        raise ValueError(
            'Target shape must contain all source shape elements'
            ' (in correct order). {} -> {}'.format(source_shape, target_shape))
    return tuple(shape)

This can be used to define unsqueeze() as a more flexible inverse of np.squeeze() compared to np.expand_dims() which can only append one singleton at a time:

def unsqueeze(
        arr,
        axis=None,
        shape=None,
        reverse=False):
    """
    Add singletons to the shape of an array to broadcast-match a given shape.

    In some sense, this function implements the inverse of `numpy.squeeze()`.


    Args:
        arr (np.ndarray): The input array.
        axis (int|Iterable|None): Axis or axes in which to operate.
            If None, a valid set axis is generated from `shape` when this is
            defined and the shape can be matched by `unsqueezing()`.
            If int or Iterable, specified how singletons are added.
            This depends on the value of `reverse`.
            If `shape` is not None, the `axis` and `shape` parameters must be
            consistent.
            Values must be in the range [-(ndim+1), ndim+1]
            At least one of `axis` and `shape` must be specified.
        shape (int|Iterable|None): The target shape.
            If None, no safety checks are performed.
            If int, this is interpreted as the number of dimensions of the
            output array.
            If Iterable, the result must be broadcastable to an array with the
            specified shape.
            If `axis` is not None, the `axis` and `shape` parameters must be
            consistent.
            At least one of `axis` and `shape` must be specified.
        reverse (bool): Interpret `axis` parameter as its complementary.
            If True, the dims of the input array are placed at the positions
            indicated by `axis`, and singletons are placed everywherelse and
            the `axis` length must be equal to the number of dimensions of the
            input array; the `shape` parameter cannot be `None`.
            If False, the singletons are added at the position(s) specified by
            `axis`.
            If `axis` is None, `reverse` has no effect.

    Returns:
        arr (np.ndarray): The reshaped array.

    Raises:
        ValueError: if the `arr` shape cannot be reshaped correctly.

    Examples:
        Let's define some input array `arr`:

        >>> arr = np.arange(2 * 3 * 4).reshape((2, 3, 4))
        >>> arr.shape
        (2, 3, 4)

        A call to `unsqueeze()` can be reversed by `np.squeeze()`:

        >>> arr_ = unsqueeze(arr, (0, 2, 4))
        >>> arr_.shape
        (1, 2, 1, 3, 1, 4)
        >>> arr = np.squeeze(arr_, (0, 2, 4))
        >>> arr.shape
        (2, 3, 4)

        The order of the axes does not matter:

        >>> arr_ = unsqueeze(arr, (0, 4, 2))
        >>> arr_.shape
        (1, 2, 1, 3, 1, 4)

        If `shape` is an int, `axis` must be consistent with it:

        >>> arr_ = unsqueeze(arr, (0, 2, 4), 6)
        >>> arr_.shape
        (1, 2, 1, 3, 1, 4)
        >>> arr_ = unsqueeze(arr, (0, 2, 4), 7)
        Traceback (most recent call last):
          ...
        ValueError: Incompatible `[0, 2, 4]` axis and `7` shape for array of\
 shape (2, 3, 4)

        It is possible to reverse the meaning to `axis` to add singletons
        everywhere except where specified (but requires `shape` to be defined
        and the length of `axis` must match the array dims):

        >>> arr_ = unsqueeze(arr, (0, 2, 4), 10, True)
        >>> arr_.shape
        (2, 1, 3, 1, 4, 1, 1, 1, 1, 1)
        >>> arr_ = unsqueeze(arr, (0, 2, 4), reverse=True)
        Traceback (most recent call last):
          ...
        ValueError: When `reverse` is True, `shape` cannot be None.
        >>> arr_ = unsqueeze(arr, (0, 2), 10, True)
        Traceback (most recent call last):
          ...
        ValueError: When `reverse` is True, the length of axis (2) must match\
 the num of dims of array (3).

        Axes values must be valid:

        >>> arr_ = unsqueeze(arr, 0)
        >>> arr_.shape
        (1, 2, 3, 4)
        >>> arr_ = unsqueeze(arr, 3)
        >>> arr_.shape
        (2, 3, 4, 1)
        >>> arr_ = unsqueeze(arr, -1)
        >>> arr_.shape
        (2, 3, 4, 1)
        >>> arr_ = unsqueeze(arr, -4)
        >>> arr_.shape
        (1, 2, 3, 4)
        >>> arr_ = unsqueeze(arr, 10)
        Traceback (most recent call last):
          ...
        ValueError: Axis (10,) out of range.

        If `shape` is specified, `axis` can be omitted (USE WITH CARE!) or its
        value is used for addiotional safety checks:

        >>> arr_ = unsqueeze(arr, shape=(2, 3, 4, 5, 6))
        >>> arr_.shape
        (2, 3, 4, 1, 1)
        >>> arr_ = unsqueeze(
        ...     arr, (3, 6, 8), (2, 5, 3, 2, 7, 2, 3, 2, 4, 5, 6), True)
        >>> arr_.shape
        (1, 1, 1, 2, 1, 1, 3, 1, 4, 1, 1)
        >>> arr_ = unsqueeze(
        ...     arr, (3, 7, 8), (2, 5, 3, 2, 7, 2, 3, 2, 4, 5, 6), True)
        Traceback (most recent call last):
          ...
        ValueError: New shape [1, 1, 1, 2, 1, 1, 1, 3, 4, 1, 1] cannot be\
 broadcasted to shape (2, 5, 3, 2, 7, 2, 3, 2, 4, 5, 6)
        >>> arr = unsqueeze(arr, shape=(2, 5, 3, 7, 2, 4, 5, 6))
        >>> arr.shape
        (2, 1, 3, 1, 1, 4, 1, 1)
        >>> arr = np.squeeze(arr)
        >>> arr.shape
        (2, 3, 4)
        >>> arr = unsqueeze(arr, shape=(5, 3, 7, 2, 4, 5, 6))
        Traceback (most recent call last):
          ...
        ValueError: Target shape must contain all source shape elements\
 (in correct order). (2, 3, 4) -> (5, 3, 7, 2, 4, 5, 6)

        The behavior is consistent with other NumPy functions and the
        `keepdims` mechanism:

        >>> axis = (0, 2, 4)
        >>> arr1 = np.arange(2 * 3 * 4 * 5 * 6).reshape((2, 3, 4, 5, 6))
        >>> arr2 = np.sum(arr1, axis, keepdims=True)
        >>> arr2.shape
        (1, 3, 1, 5, 1)
        >>> arr3 = np.sum(arr1, axis)
        >>> arr3.shape
        (3, 5)
        >>> arr3 = unsqueeze(arr3, axis)
        >>> arr3.shape
        (1, 3, 1, 5, 1)
        >>> np.all(arr2 == arr3)
        True
    """
    # calculate `new_shape`
    if axis is None and shape is None:
        raise ValueError(
            'At least one of `axis` and `shape` parameters must be specified.')
    elif axis is None and shape is not None:
        new_shape = unsqueezing(arr.shape, shape)
    elif axis is not None:
        if isinstance(axis, int):
            axis = (axis,)
        # calculate the dim of the result
        if shape is not None:
            if isinstance(shape, int):
                ndim = shape
            else:  # shape is a sequence
                ndim = len(shape)
        elif not reverse:
            ndim = len(axis) + arr.ndim
        else:
            raise ValueError('When `reverse` is True, `shape` cannot be None.')
        # check that axis is properly constructed
        if any([ax < -ndim - 1 or ax > ndim + 1 for ax in axis]):
            raise ValueError('Axis {} out of range.'.format(axis))
        # normalize axis using `ndim`
        axis = sorted([ax % ndim for ax in axis])
        # manage reverse mode
        if reverse:
            if len(axis) == arr.ndim:
                axis = [i for i in range(ndim) if i not in axis]
            else:
                raise ValueError(
                    'When `reverse` is True, the length of axis ({})'
                    ' must match the num of dims of array ({}).'.format(
                        len(axis), arr.ndim))
        elif len(axis) + arr.ndim != ndim:
            raise ValueError(
                'Incompatible `{}` axis and `{}` shape'
                ' for array of shape {}'.format(axis, shape, arr.shape))
        # generate the new shape from axis, ndim and shape
        new_shape = []
        i, j = 0, 0
        for l in range(ndim):
            if i < len(axis) and l == axis[i] or j >= arr.ndim:
                new_shape.append(1)
                i += 1
            else:
                new_shape.append(arr.shape[j])
                j += 1

    # check that `new_shape` is consistent with `shape`
    if shape is not None:
        if isinstance(shape, int):
            if len(new_shape) != ndim:
                raise ValueError(
                    'Length of new shape {} does not match '
                    'expected length ({}).'.format(len(new_shape), ndim))
        else:
            if not all([new_dim == 1 or new_dim == dim
                        for new_dim, dim in zip(new_shape, shape)]):
                raise ValueError(
                    'New shape {} cannot be broadcasted to shape {}'.format(
                        new_shape, shape))

    return arr.reshape(new_shape)

Using these, one can write:

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))

arr3 = unsqueeze(arr2, (0, 2, 4))

arr1 + arr3
# now this works because it has the correct shape

arr3 = unsqueeze(arr2, shape=arr1.shape)
arr1 + arr3
# this also works because the shape can be expanded unambiguously

So dynamic broadcast can now happen, and this is consistent with the behavior of keepdims:

import numpy as np

axis = (0, 2, 4)
arr1 = np.arange(2 * 3 * 4 * 5 * 6).reshape((2, 3, 4, 5, 6))
arr2 = np.sum(arr1, axis, keepdims=True)
arr3 = np.sum(arr1, axis)
arr3 = unsqueeze(arr3, axis)
np.all(arr2 == arr3)
# : True

Effectively, this extends np.expand_dims() to handle more complex scenarios.

Improvements over this code are obviously more than welcome.

norok2
  • 25,683
  • 4
  • 73
  • 99
  • It's up to you to prove that the logic is correct and unambiguous, at least for your own use cases. – hpaulj Sep 27 '17 at 08:14
  • I have extended `unsqueeze()` to handle an `axis` parameters to mimic `np.expand_dims()` for when its a sequence of ints, and kept the *automagic* expansion to when `axis` is not specified, additionally warning users when *automagic* results are not uniquely defined. – norok2 Sep 27 '17 at 12:12