Vector-Scalar Example

The simplest example I can think of is the subtraction of a scalar from a vector. This operation is commonly used to center, or shift, a vector with respect to a given value. Let’s say we have the vector \(x\) of a given dimension \(d\) and want to subtract the value \(1\) from each element. Assuming we have the values of \(x\) in the numpy array x, we can use a for-loop and do:

y = x.copy()
for i in range(len(x)):
  y[i] = x[i] - 1

or, we can do:

y = x - 1

The last piece of code is shorter and more intuitive since that’s how we also write that operation in paper; but how does numpy know how to subtract 1 from each element in x? That’s where broadcasting comes in; the image below shows how it happens in this case,

The value \(1\) was repeated, or broadcasted, along the longest dimension of \(x\) and then subtracted from each element of \(x\). Broadcasting happened this way because of the shapes, or dimensions, of \(1\) and \(x\); had their shapes been different, we would’ve had to manipulate these in order to perform the subtraction. In fact, the trick to understand broadcasting lays upon understanding the shapes of the terms we’re dealing with, and how to modify them.

Matrix-Vector Example

Let’s now increase the dimension of the objects and subtract a vector \(\mu\) from a matrix \(X\) and call the difference \(dX\); this operation is equivalent to centering the data in a matrix around a point. Each row of \(X\) is composed of a d-dimensional vector \(x_i\), and the centering vector \(\mu\) shares the same dimension of \(x_i\). In matrix notation, the data set is

\[ X = \left[\begin{array}{c}x_1^T\\ x_2^T \\ \vdots \\ x_n^T \end{array}\right] \]

and it’s centered as

\[ dX = X - 1_{n} \mu^T = \left[\begin{array}{c}x_1^T - \mu^T\\ x_2^T - \mu^T \\ \vdots \\ x_n^T - \mu^T \end{array} \right] \]

The term \(1_{n}\) is a vector with \(n\) values of \(1\). Now, assuming that we have numpy arrays with the matrix X and vector mu, we obtain the centered matrix dX with the line

dX = X - mu

The way the code is written is very similar to how we wrote the operation in the equation above; this similarity is one of the beautiful and powerful features of numpy broadcasting. The above computation happens assuming that the shape of X is \((n, d)\), and mu is a 1D array of shape \((d, )\), or a 2D array of shape \((1, d)\). Why can we use a 1D or 2D array for mu? This is because in numpy 1D arrays are considered row vectors when using them in operations with higher dimensional arrays. Furthermore, the rightmost value of the shape of an array corresponds to the number of columns of the array. In this case the matrix \(X\) and vector \(\mu\) match the number of elements \(d\) in the column dimension. Broadcasting repeats the vector \(\mu\) \(n\) times along the row-dimension. Here’s a sketch of that:

Matrix-Matrix Example

Now, let’s say you want to center the matrix \(X\) around \(m\) different vectors \(\mu_j\) contained in a matrix \(M\) that looks like this

\[ M = \left[\begin{array}{c}\mu_1^T\\ \mu_2^T \\ \vdots \\ \mu_m^T \end{array} \right] \]

This type of centering comes up in algorithms like K-Means clustering or Expectation-Maximization. In this case we can use broadcasting by arranging \(M\) in such a way that \(X\) is repeated along a dimension of \(M\) from which \(\mu_j\) is subtracted. We can visualize this in the following sketch,

Now, how do we arrange the shape of \(M\) in the form shown above? That’s where array shapes and reshaping come into play.

Array Shapes

I already mentioned that the first number on the right of an array shape corresponds to the number of columns in the array; now, the second number from the right corresponds to the rows of the array, and the third number corresponds to its number of layers. If we consider \(M\) to be in a 3-dimensional space, counting from the right, the first number corresponds to the x-axis, the second to the y-axis and the third to the z-axis. Knowing this, we can add extra dimensions to \(M\) to shape it in a way that’s useful for us. If its shape is \((m, d)\), it means that we have \(m\) rows and \(d\) columns, but if we add another dimension and make the shape \((m, 1, d)\), we now have \(m\) layers, \(1\) row and \(d\) columns; similarly, a shape of \((m, d, 1)\) provides us with \(m\) layers, \(d\) rows and \(1\) column. These shape transformations look like this:

You increase the dimension of an array by slicing with np.newaxis or None in the position of the dimension you want to add. Assuming that M is a numpy array with shape \((m, d)\), you obtain the shapes

• \((1, m, d)\) with: M[None, :, :]
• \((m, 1, d)\) with: M[:, None, :]
• \((m, d, 1)\) with: M[:, :, None]

Broadcasting in multiple dimensions

From the figure above, the second item in the previous list corresponds to the reshaped version of M we need. Thus, the broadcasted version of \(X - \mu_j\) would be computed as

dX = X - M[:, None, :]

In this operation, numpy expands the shape of X to \((1, n, d)\), and broadcasts its values \(m\) times along the first dimension producing an \((m, n, d)\) array; then it repeats \(n\) times the reshaped version of M (with shape \((m, 1, d)\)) along its second dimension, producing another \((m, n, d)\) array; finally, it performs element-wise subtraction to create the array dX with shape \((m, n, d)\).

If we are interested in the distance between each point \(x_i\) and vector \(\mu_j\) we just calculate the 2-norm of the array dX along its column dimension; we do this because the difference \(x_i - \mu_j\) is stored in each row of dX, and we’d have to move across its columns in order to sum the squared values of the elements of \(x_i - \mu_j\). The code that does this is:

distance = np.linalg.norm(dX, axis=2)

We didn’t specify the 2-norm in the previous function, because this norm is a default parameter. For a different norm, we would have to include the parameter ord; for example, ord=1 or ord=inf.

Since we calculated the norm across the column dimension of dX, this dimension will vanish and the new array will keep the shape of the remaining dimensions, this is \((m, n)\); therefore, the array distance will have \(m\) rows and \(n\) columns. The resulting distance between the vectors \(x_i\) and \(\mu_j\) is stored in the ijth position of the distance array.