Broadcasting Explained

1 Introduction

In the context of operations involving 2 (or more) arrays, “broadcasting” refers to recycling array dimensions without allocating additional memory or making needles copies. This is considerably faster and more memory-efficient than ’s regular dimensions repetition (and similar) mechanisms.

This page explains the concept of “broadcasting” in more detail. A good understanding of atomic and recursive arrays in base is somewhat essential to follow this page.

2 What is broadcasting and why is it needed?

2.1 Example case

Let’s start with a simple example.
Consider the matrices x and y:

x <- array(1:20, c(4, 5))
y <- array(1:5*10, c(1, 5))
print(x)
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    1    5    9   13   17
#> [2,]    2    6   10   14   18
#> [3,]    3    7   11   15   19
#> [4,]    4    8   12   16   20
print(y)
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   10   20   30   40   50

Suppose one wishes to compute the element-wise addition of these 2 arrays.

This won’t work in base :

x + y
Error in x + y : non-conformable arrays

When computing the element-wise sum of these arrays, one or both of them need to be recycled so that they are equal size, in order to compute the element-wise computation.

In this case, matrix y needs its single row to be recycled 4 times, making y the same size as x, and thus conformable.

provides linear vector recycling, but not recycling of array dimensions. Instead, in base we need to replicate (and thus copy) the array dimensions. This can be done manually, like so:

x + y[rep(1L, 4L),]
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   11   25   39   53   67
#> [2,]   12   26   40   54   68
#> [3,]   13   27   41   55   69
#> [4,]   14   28   42   56   70

2.2 What is the problem?

There are several problems.

First and foremost, when the arrays become larger, the replicated arrays become larger, and if the arrays become too large, you may require more memory than is available in your current system, resulting in a message like the following:

> Error: cannot allocate vector of size

The problem isn’t limited to available memory. As the required memory to allocate increase, the speed decreases. And let’s not forget that computational inefficiency in general is also bad for the environment.

Additionally, a solution like x + y[rep(1L, 4L),] is not easily scalable for other arrays when the dimensions of x and y are not known a-priori.

2.3 Introducing Broadcasting

In an operation like x + y[rep(1L, 4L),], y is replicated to become the same size as x. The thing is, physical replication - and thus copying - of the dimensions of an array should not be necessary; arrays only need to be recycled virtually.

Virtual recycling does not actually physically replicate arrays. Instead, nested loops in ‘C’ and ‘C++’ are used to simulate a recycled array.
This is similar to how recycles regular (i.e. dimensionless) vectors.
Virtual recycling requires no additional memory (apart from allocating the final end result); it is much faster and much more memory efficient than replicating dimensions.

And that is what broadcasting does: broadcasting provides fast virtual recycling of array dimensions in the context of element-wise operations involving 2 (or more) arrays without allocating additional memory. Broadcasting in this package is also scalable to arrays of any dimensions (up to 16 dimensions).

In the earlier example, we used:

x + y[rep(1L, 4L),]
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   11   25   39   53   67
#> [2,]   12   26   40   54   68
#> [3,]   13   27   41   55   69
#> [4,]   14   28   42   56   70

To compute the element-wise addition using broadcasting through the ‘broadcast’ package, we can do the following:

library(broadcast)

bc.num(x, y, "+")
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   11   25   39   53   67
#> [2,]   12   26   40   54   68
#> [3,]   13   27   41   55   69
#> [4,]   14   28   42   56   70

The result is the same (as it should). But as the size of the resulting array increases, the broadcasted functions become more and more efficient in terms of both speed and memory, in comparison with base approaches.

Technically, one can also use t(t(x) + drop(y)), but that has a similar problem: one needs to perform transposition twice, which makes unnecessary copies, is slow, and is - again - not scalable to arrays where the dimensions are not known a-priori.

Benchmarks can be found on the website.

3 Illustrating Broadcasting

3.1 Retracing the first example

In the previous example the following arrays were used:

x <- array(1:20, c(4, 5))
y <- array(1:5*10, c(1, 5))

Let’s see what happens if we sum them together using broadcasting:

bc.num(x, y, "+")
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   11   25   39   53   67
#> [2,]   12   26   40   54   68
#> [3,]   13   27   41   55   69
#> [4,]   14   28   42   56   70

What happens here is that x remains the same, but row vector y is virtually recycled 4 times, but without requiring 4 times more memory. This is also illustrated here:

x[1, 1] + y[1, 1]	x[1, 2] + y[1, 2]	x[1, 3] + y[1, 3]	x[1, 4] + y[1, 4]	x[1, 5] + y[1, 5]
x[2, 1] + y[1, 1]	x[2, 2] + y[1, 2]	x[2, 3] + y[1, 3]	x[2, 4] + y[1, 4]	x[2, 5] + y[1, 5]
x[3, 1] + y[1, 1]	x[3, 2] + y[1, 2]	x[3, 3] + y[1, 3]	x[3, 4] + y[1, 4]	x[3, 5] + y[1, 5]
x[4, 1] + y[1, 1]	x[4, 2] + y[1, 2]	x[4, 3] + y[1, 3]	x[4, 4] + y[1, 4]	x[4, 5] + y[1, 5]

3.2 Perpendicular vectors

Let’s now consider another scenario. We now take again 2 matrices x and y, but this time x is a column vector (i.e. a matrix with 1 column and multiple rows), and y is a row vector (i.e. a matrix with 1 row and multiple columns):

x <- array(1:5, c(5, 1))
y <- array(1:5 * 10, c(1, 5))
print(x)
#>      [,1]
#> [1,]    1
#> [2,]    2
#> [3,]    3
#> [4,]    4
#> [5,]    5
print(y)
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   10   20   30   40   50

Computing the broadcasted element-wise sum of x and y produces the following:

bc.num(x, y, "+")
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   11   21   31   41   51
#> [2,]   12   22   32   42   52
#> [3,]   13   23   33   43   53
#> [4,]   14   24   34   44   54
#> [5,]   15   25   35   45   55

So what exactly does broadcasting compute here?

The following table illustrates what broadcasting does, conceptually:

x[1, 1] + y[1, 1]	x[1, 1] + y[1, 2]	x[1, 1] + y[1, 3]	x[1, 1] + y[1, 4]	x[1, 1] + y[1, 5]
x[2, 1] + y[1, 1]	x[2, 1] + y[1, 2]	x[2, 1] + y[1, 3]	x[2, 1] + y[1, 4]	x[2, 1] + y[1, 5]
x[3, 1] + y[1, 1]	x[3, 1] + y[1, 2]	x[3, 1] + y[1, 3]	x[3, 1] + y[1, 4]	x[3, 1] + y[1, 5]
x[4, 1] + y[1, 1]	x[4, 1] + y[1, 2]	x[4, 1] + y[1, 3]	x[4, 1] + y[1, 4]	x[4, 1] + y[1, 5]
x[5, 1] + y[1, 1]	x[5, 1] + y[1, 2]	x[5, 1] + y[1, 3]	x[5, 1] + y[1, 4]	x[5, 1] + y[1, 5]

4 Broadcasting limitations and rules

The examples shown so far used matrices. But ‘broadcast’ is not restricted to only handling matrices; this package can handle any array, provided it does not have more than 16 dimensions. Moreover, any number of dimensions (provided <= 16) can be broadcasted simultaneously.

The examples show so far only used addition. But ‘broadcast’ supports a wide variety of binary operations, inlcuding relational-, arithmetic-, Boolean-, and string operations.

4.1 Normalization

Before broadcasting occurs, the number of dimensions of x and y are normalized, such that ndim(x) and ndim(y) are the same ( the ndim(x) function is the same as length(dim(x))).

This normalization achieved by appending dimensions of size 1 to the array with the smaller number of dimensions until both array have the same number of dimensions.

For example, if dim(x) = c(3, 4, 3) and dim(y) = c(3, 4), then internally dim(y) is changed to c(3, 4, 1).

4.2 Conformability

Broadcasted operations only work when 2 arrays are conformable for broadcasting.

Consider again 2 arrays x and y, and their dimensions given by dim(x) and dim(y), respectively.
‘broadcast’ goes through dim(x) and dim(y), from left (i.e. first dimension / rows) to right (i.e. the last dimension), and checks for each axis i if at least one of the following conditions is TRUE:

dim(x)[i] and dim(y)[i] are equal;
either dim(x)[i] or dim(y)[i] is 1;
either dim(x)[i] or dim(y)[i] is non-existing (see also the “Normalization” sub-section above)

if at least one of the above conditions is true for each and every one of the dimensions of x and y, the 2 arrays are compatible for broadcasted operations. If not, they are not compatible, and attempting to perform a broadcasted operations results in an error.

To illustrate, let’s check 2 arrays to see if they are compatible:

x <- array(rnorm(10), c(10, 1, 9, 6))
y <- array(rnorm(10), c(10, 5, 1))

i	dim(x)[i]	dim(y)[i]	compatible	reason
1	10	10	TRUE	equal
2	1	5	TRUE	either is 1
3	9	1	TRUE	either is 1
4	6	NA	TRUE	either is missing

As shown in the data.frame above, all dimensions of x and y are compatible.
Therefore, broadcasted binary operations involving x and y can be performed.

Now let’s look at another 2 arrays, and see if they are compatible:

x <- array(rnorm(10), c(10, 1, 9))
y <- array(rnorm(10), c(10, 5, 2, 6))

i	dim(x)[i]	dim(y)[i]	compatible	reason
1	10	10	TRUE	equal
2	1	5	TRUE	either is 1
3	9	2	FALSE	not equal & neither is 1/missing
4	NA	6	TRUE	either is missing

One of the dimensions, namely the third dimension, is not compatible. Therefore, these 2 arrays are not conformable for broadcasted operations.

4.3 Which dimensions are broadcasted

In the context of an operation involving exactly 2 arrays, this sub-section explains which dimensions of which array is broadcasted, and which dimension of which array is left as-is.
There are 3 scenarios relevant for this explanation.

Scenario 1:
For some dimension i, dim(x)[i] and dim(y)[i] are equal.
In this scenario, no broadcasting needs to occur on that dimension for either array.

Scenario 2:
For some dimension i, suppose dim(x)[i] == 1 and dim(y)[i] > 1.
Then dimension dim(x)[i] is broadcasted to size dim(y)[i].
Similarly, if dim(x)[i] > 1 and dim(y)[i] == 1, dimension dim(y)[i] is broadcasted to size dim(x)[i].

Scenario 3:
Finally, suppose for some dimension i, dim(x)[i] > 1, and dim(y)[i] is non-existing.
I.e., x has more dimensions than y.
In this case, the missing dimension in y is replaced with a dimensions of size 1, and scenario 2 is used.
Similarly, if dim(y)[i] > 1 and dim(x)[i] is non-existing, the missing dimension in x is replaced with a dimension of size 1 and we get back to scenario 2.

It does not matter how much difference there is between the number of dimensions of x and the number of dimensions of y: any number of missing dimensions will be replaced with 1, as long as the total number of dimensions for each array does not exceed 16.

Illustration

To illustrate, let’s check 2 conformable arrays, and see which dimension of which array is broadcasted by how much

x <- array(rnorm(10), c(10, 1, 1, 9, 6))
y <- array(rnorm(10), c(10, 1, 5, 1))

i	dim(x)[i]	dim(y)[i]	broadcasted	reason	note
1	10	10	neither	dims equal
2	1	1	neither	dims equal
3	1	5	x	dim(x)[i] == 1	dim(x)[i] recycled to size 5
4	9	1	y	dim(y)[i] == 1	dim(y)[i] recycled to size 9
5	6	NA	y	dim(y)[i] is missing	dim(y) extended with a 1; dim(y)[i] recycled to size 6

5 Orthogonal Arrays

In the documentation of ‘broadcast’, the reader may come across the term “orthogonal arrays” here and there.
“orthogonal” is a term that can mean a great many things; this section therefore explains what “orthogonal” means in the context of broadcasting.

Consider a column-vector x and row-vector y:

x <- array(1:5, c(5, 1))
y <- array(1:5*10, c(1, 5))
print(x)
#>      [,1]
#> [1,]    1
#> [2,]    2
#> [3,]    3
#> [4,]    4
#> [5,]    5
print(y)
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   10   20   30   40   50

These vectors are (in a very literal sense) perpendicular to each other, and thus in that sense orthogonal to each other.
This type of orthogonality can be extended to arrays as follows.
For the purposes of broadcasting, any pair of arrays, x and y, can be referred to as being “orthogonal” if all of the following holds for that pair of arrays:

For every dimension index i, dim(x)[i] != dim(y)[i]
For every dimensions index i,dim(x)[i] is 1 or missing, OR dim(y)[i] is 1 or missing.

5.1 Illustration

To illustrate, here are 2 arrays that are fully orthogonal:

x <- array(1:10, c(10, 1, 8, 1, 10))
y <- array(1:10, c(1, 9, 1, 9, 1, 11))

i	dim(x)[i]	dim(y)[i]	orthogonal	reason
1	10	1	TRUE	dim(x)[i] != dim(y)[i]
2	1	9	TRUE	dim(x)[i] != dim(y)[i]
3	8	1	TRUE	dim(x)[i] != dim(y)[i]
4	1	9	TRUE	dim(x)[i] != dim(y)[i]
5	10	1	TRUE	dim(x)[i] != dim(y)[i]
6	NA	11	TRUE	one of them missing

And the following 2 arrays are not orthogonal:

x <- array(1:10, c(10, 1, 8, 1, 10))
y <- array(1:10, c(10, 9, 1, 9, 10))

i	dim(x)[i]	dim(y)[i]	orthogonal	reason
1	10	10	FALSE	dim(x)[i] == dim(y)[i]
2	1	9	TRUE	dim(x)[i] != dim(y)[i]
3	8	1	TRUE	dim(x)[i] != dim(y)[i]
4	1	9	TRUE	dim(x)[i] != dim(y)[i]
5	10	10	FALSE	dim(x)[i] == dim(y)[i]

The above arrays are not orthogonal, due to dimensions 1 and 5.

6 Broadcasting vs outer

The base outer() function performs outer computations on vectors and arrays.

Earlier in this guide, we used the following computation:

x <- array(1:5, c(5, 1))
y <- array(1:5*10, c(1, 5))

bc.num(x, y, "+")
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   11   21   31   41   51
#> [2,]   12   22   32   42   52
#> [3,]   13   23   33   43   53
#> [4,]   14   24   34   44   54
#> [5,]   15   25   35   45   55

This is equivalent to using the outer() function as follows:

outer(as.vector(x), as.vector(y), "+")
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   11   21   31   41   51
#> [2,]   12   22   32   42   52
#> [3,]   13   23   33   43   53
#> [4,]   14   24   34   44   54
#> [5,]   15   25   35   45   55

As you can see, performing a broadcasted operation with a row-vector and a column-vector (in that order) is equivalent to using outer() on 2 dimensionless vectors.

The outer() approach has the same problems as the x[, rep(1L, 5L)] + y[rep(1L, 5L),] approach: it consumes a lot of unnecessary memory, and is slow.

Broadcasting can thus be used as a more efficient alternative to outer() in case of vector inputs, simply by making the inputs orthogonal. But please don’t confuse outer computations with broadcasting: broadcasting is efficient array recycling for element-wise computations.