Take 25% off Julia for Data Analysis by entering fcckaminski into the discount code box at checkout at manning.com.

We will discuss how broadcasting works by using Anscombe’s quartet data. However, let us start with explaining of broadcasting on some toy examples.

Syntax and meaning of broadcasting in Julia

An important design rule of the Julia language is that it uses the approach that definitions of functions follow the rules of mathematics. The multiplication operator * uses matrix multiplication rules. Therefore, the following code follows matrix multiplication rules:

 
 julia> x = [1 2 3]
 1×3 Matrix{Int64}:
  1  2  3
  
 julia> y = [1, 2, 3]
 3-element Vector{Int64}:
  1
  2
  3
  
 julia> x * y
 1-element Vector{Int64}:
  14
  

The operation works as we multiply x which is bound to a 1×3 matrix and a 3-element vector y and in Julia vectors are always interpreted as columnar.

You might ask then how one should multiply two vectors elementwise, an operation known as Hadamard product in mathematics. Clearly this is not possible just with the * operator as it does standard matrix multiplication:

 
 julia> a = [1, 2, 3]
 3-element Vector{Int64}:
  1
  2
  3
  
 julia> b = [4, 5, 6]
 3-element Vector{Int64}:
  4
  5
  6
  
 julia> a * b
 ERROR: MethodError: no method matching *(::Vector{Int64}, ::Vector{Int64})
  

You get an error as multiplication of a vector by a vector is not a valid mathematical operation. What we need to do is to broadcast the multiplication. In Julia adding broadcasting to an operator is easy. You just prefix it with a dot (.) like this:

 
 julia> a .* b
 3-element Vector{Int64}:
   4
  10
  18
  

When a broadcasted operator like .* is used Julia iterates elements of passed collections, in our case vectors a and b, and applies the operator after the dot, in our case *, elementwise. Therefore, in this case broadcasting result is the same as produced by the following operations:

 
 julia> map(*, a, b)
 3-element Vector{Int64}:
   4
  10
  18
  
 julia> [a[i] * b[i] for i in eachindex(a, b)]
 3-element Vector{Int64}:
   4
  10
  18
  

In the map example we are passing two collections (instead of only one as we did before when we explained how map works) and the passed function, * in this case, is applied iteratively elementwise to them until one of the collections gets exhausted.

In the comprehension example it is worth to comment on the eachindex(a, b) expression which produces:

 
 julia> eachindex(a, b)
 Base.OneTo(3)
  

The eachindex function produces indices that could be used to index into both a and b arguments passed to it. In this case these are just integers from 1 to 3. This means that you can index both a and b vectors with these values; for example, the following indexing expressions are valid: a[1], b[2], a[3], but a[0] or b[4] would not be valid as they are not in the range specified by Base.OneTo(3).

If sizes of a and b would not match, we get an error:

 
 julia> eachindex([1, 2, 3], [4, 5])
 ERROR: DimensionMismatch("all inputs to eachindex must have the same indices, got Base.OneTo(3) and Base.OneTo(2)")
  

This is a very important difference from the map function which does not use the eachindex function internally, but instead iterates collections until either of them is exhausted as we have explained above:

 
 julia> map(*, [1, 2, 3], [4, 5])
 2-element Vector{Int64}:
   4
  10
  

Using broadcasting just like the eachindex function checks if the dimensions of the passed objects match:

 
 julia> [1, 2, 3] .* [4, 5]
 ERROR: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 3 and 2")
  

Expansion of length-1 dimensions in broadcasting

There is one exception to the rule that dimensions of all collections taking part in broadcasting must match. This additional rule states that single element dimensions get expanded to match the size of the other collection:

 
 julia> [1, 2, 3] .^ [2]
 3-element Vector{Int64}:
   1
   4
   9
  

You might ask why single element dimension gets expanded? The reason is practical: in most cases when your collection has a single element in some dimension you want it to get expanded. Here is a first example:

 
 julia> [1, 2, 3] .^ 2
 3-element Vector{Int64}:
  1
  4
  9
  

Here we are calculating a square of elements of a vector. Since 2 is a scalar, it is interpreted as having size 1 in each dimension. Most people agree that in such a case dimension expansion should happen. You will see the same behavior both in Python and R.

Now let us consider a second example:

 
 julia> [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] .* [1 2 3 4 5 6 7 8 9 10]
 10×10 Matrix{Int64}:
   1   2   3   4   5   6   7   8   9   10
   2   4   6   8  10  12  14  16  18   20
   3   6   9  12  15  18  21  24  27   30
   4   8  12  16  20  24  28  32  36   40
   5  10  15  20  25  30  35  40  45   50
   6  12  18  24  30  36  42  48  54   60
   7  14  21  28  35  42  49  56  63   70
   8  16  24  32  40  48  56  64  72   80
   9  18  27  36  45  54  63  72  81   90
  10  20  30  40  50  60  70  80  90  100
  

Here we have created a multiplication table. The reason why the specified operation worked is that we had a 10-element vector having 1 column and 10 rows, and a 10-element matrix having 1 row and 10 columns. In this case dimension expansion happened both for the left- and the right-hand side of the operation.

The technique we have shown above is used very often in practice to get a Cartesian product of all inputs. For instance, in part 2 of the book you will learn that in DataFrames.jl when you write "x" => sum you ask the package to apply the sum function to the column "x" of the data frame. A common task is that we want to apply several functions to several columns of a data frame. Using broadcasting this can be written concisely as:

 
 julia> ["x", "y", "z"] .=> [sum minimum maximum]
 3×3 Matrix{Pair{String, _A} where _A}:
  "x"=>sum  "x"=>minimum  "x"=>maximum
  "y"=>sum  "y"=>minimum  "y"=>maximum
  "z"=>sum  "z"=>minimum  "z"=>maximum
  

Such an expression would ask for computation of sum, minimum, and maximum for columns "x", "y", and "z".

You now know that you can add a dot (.) before any operator to broadcast it. How about functions that are not operators? Here you also use a dot (.), but this time you suffix it after the function name. Here is an example:

 
 julia> abs.([1, -2, 3, -4])
 4-element Vector{Int64}:
  1
  2
  3
  4
  

Let me stress that just applying the abs function to a vector results in an error:

 
 julia> abs([1, 2, 3])
 ERROR: MethodError: no method matching abs(::Vector{Int64})
  

The reason is the same as before: absolute value is mathematically defined for numbers but is not defined for vectors. Of course, you can conveniently also apply broadcasting to a function taking multiple arguments. For example, the string function concatenates its arguments into a single string:

 
 julia> string(1, 2, 3)
 "123"
  

If we use broadcasting on this function, we get for instance the following result:

 
 julia> string.("x", 1:10)
 10-element Vector{String}:
  "x1"
  "x2"
  "x3"
  "x4"
  "x5"
  "x6"
  "x7"
  "x8"
  "x9"
  "x10"
  

Here we expanded the dimension of the scalar "x" to match the length of the 1:10 range. Such an operation is quite common when one wants to automatically generate names for some objects, for example, columns of some data frame or file names in some folder.

It is important to highlight here that prefixing a dot before the operator or suffixing it to a function name is a fully general solution. It is not a hardcoded functionality of some specific predefined operations. You can use broadcasting with any custom function. For instance:

 
 julia> f(i::Int) = string("got integer ", i)
 f (generic function with 2 methods)
  
 julia> f(s::String) = string("got string ", s)
 f (generic function with 2 methods)
  
 julia> f.([1, "1"])
 2-element Vector{String}:
  "got integer 1"
  "got string 1"
  

Here we have defined two methods for the f function. As you can see by writing f., we have automatically broadcasted it without having to define anything extra.

Protection of collections from being broadcasted over

Before we go back to our Anscombe’s quartet data let me comment on one very common case. What should one do if one does not want to broadcast some collection but to force its reuse along all dimensions as if it were a scalar? To understand this issue let me introduce the in function first.

 
 julia> in(1, [1, 2, 3])
 true
  
 julia> in(4, [1, 2, 3])
 false
  

 
 julia> 1 in [1, 2, 3]
 true
  
 julia> 4 in [1, 2, 3]
 false
  

Now imagine I have a long vector of values I want to check if they are contained in some vector. When you try the test without broadcasting it does not work like you probably expected:

 
 julia> in([1, 3, 5, 7, 9], [1, 2, 3, 4])
 false
  

The problem is that the vector [1, 3, 5, 7, 9] is not an element of vector [1, 2, 3, 4] so we get false. For a reference let us test the scenario when we put the [1, 3, 5, 7, 9] vector into the collection in which we look for it:

 
 julia> in([1, 3, 5, 7, 9], [1, 2, 3, 4, [1, 3, 5, 7, 9]])
 true
  

As expected, this time the in test returned true.

Going back to the original test note that broadcasting does not seem to work either:

 
 julia> in.([1, 3, 5, 7, 9], [1, 2, 3, 4])
 ERROR: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 5 and 4")
  

How should we resolve this issue? The solution is to wrap the vector that we want to be expanded with Ref. In this way we will protect this object from being iterated over. Instead, it will be unwrapped from Ref and treated by broadcasting as if it were a scalar:

 
 julia> in.([1, 3, 5, 7, 9], Ref([1, 2, 3, 4]))
 5-element BitVector:
  1
  1
  0
  0
  0
  

This time we got the expected result.

In the output note that Boolean true is printed as 1 and Boolean false is printed as 0. This choice of display was made to allow for more convenient visual inspection of large matrices containing Boolean values. To see why this is useful consider if we wanted to check which entries in a multiplication table are odd using the isodd function:

 
 julia> isodd.([1, 2, 3, 4, 5, 6, 7, 8, 9, 10] .* [1 2 3 4 5 6 7 8 9 10])
 10×10 BitMatrix:
  1  0  1  0  1  0  1  0  1  0
  0  0  0  0  0  0  0  0  0  0
  1  0  1  0  1  0  1  0  1  0
  0  0  0  0  0  0  0  0  0  0
  1  0  1  0  1  0  1  0  1  0
  0  0  0  0  0  0  0  0  0  0
  1  0  1  0  1  0  1  0  1  0
  0  0  0  0  0  0  0  0  0  0
  1  0  1  0  1  0  1  0  1  0
  0  0  0  0  0  0  0  0  0  0
  

In this example I have additionally shown you that broadcasting operations can be chained together in a single expression. In this case we broadcasted both the multiplication * and the isodd function.

For a reference let me show you how this matrix would be displayed if we changed its element type to Any:

 
 julia> Matrix{Any}(isodd.([1, 2, 3, 4, 5, 6, 7, 8, 9, 10] .*
                           [1 2 3 4 5 6 7 8 9 10]))
 10×10 Matrix{Any}:
   true  false   true  false   true  false   true  false   true  false
  false  false  false  false  false  false  false  false  false  false
   true  false   true  false   true  false   true  false   true  false
  false  false  false  false  false  false  false  false  false  false
   true  false   true  false   true  false   true  false   true  false
  false  false  false  false  false  false  false  false  false  false
   true  false   true  false   true  false   true  false   true  false
  false  false  false  false  false  false  false  false  false  false
   true  false   true  false   true  false   true  false   true  false
  false  false  false  false  false  false  false  false  false  false
  

This time true and false is printed to avoid potential confusion with integer 1 and 0 that could be potentially stored in this matrix. In my opinion, however, analyzing such a printout is less convenient than before.

To practice what we have learned try the following exercise that is a quite common task when processing data.

Analyzing Anscombe’s quartet data using broadcasting

Now we are ready to load up our Anscombe’s quartet data. Let us initialize the aq variable first:

 
 julia> aq = [10.0   8.04  10.0  9.14  10.0   7.46   8.0   6.58
               8.0   6.95   8.0  8.14   8.0   6.77   8.0   5.76
              13.0   7.58  13.0  8.74  13.0  12.74   8.0   7.71
               9.0   8.81   9.0  8.77   9.0   7.11   8.0   8.84
              11.0   8.33  11.0  9.26  11.0   7.81   8.0   8.47
              14.0   9.96  14.0  8.1   14.0   8.84   8.0   7.04
               6.0   7.24   6.0  6.13   6.0   6.08   8.0   5.25
               4.0   4.26   4.0  3.1    4.0   5.39  19.0  12.50
              12.0  10.84  12.0  9.13  12.0   8.15   8.0   5.56
               7.0   4.82   7.0  7.26   7.0   6.42   8.0   7.91
               5.0   5.68   5.0  4.74   5.0   5.73   8.0   6.89]
 11×8 Matrix{Float64}:
  10.0   8.04  10.0  9.14  10.0   7.46   8.0   6.58
   8.0   6.95   8.0  8.14   8.0   6.77   8.0   5.76
  13.0   7.58  13.0  8.74  13.0  12.74   8.0   7.71
   9.0   8.81   9.0  8.77   9.0   7.11   8.0   8.84
  11.0   8.33  11.0  9.26  11.0   7.81   8.0   8.47
  14.0   9.96  14.0  8.1   14.0   8.84   8.0   7.04
   6.0   7.24   6.0  6.13   6.0   6.08   8.0   5.25
   4.0   4.26   4.0  3.1    4.0   5.39  19.0  12.5
  12.0  10.84  12.0  9.13  12.0   8.15   8.0   5.56
   7.0   4.82   7.0  7.26   7.0   6.42   8.0   7.91
   5.0   5.68   5.0  4.74   5.0   5.73   8.0   6.89
  
 julia> using Statistics
  

We will perform two tasks using broadcasting: calculation of the mean of every variable and calculation of the coefficient of determination.

We first start with calculating the mean of the columns of the aq matrix. When thinking how to do it we notice that we need to broadcast the mean function over a collection of columns of aq. Fortunately we know that the eachcol function gives us such a collection, therefore we can write:

 
 julia> mean.(eachcol(aq))
 8-element Vector{Float64}:
  9.0
  7.500909090909093
  9.0
  7.500909090909091
  9.0
  7.500000000000001
  9.0
  7.50090909090909
  

Note the dot (.) after mean. If we would forget to write it, we would get:

 
 julia> mean(eachcol(aq))
 11-element Vector{Float64}:
   8.6525
   7.4525
  10.47125
   8.56625
   9.35875
  10.492500000000001
   6.3375
   7.03125
   9.71
   6.92625
   5.755000000000001
  

The result is just a mean of 8 vectors constituting columns of aq matrix, clearly not what we wanted.

As a second application let us rewrite the function calculating the coefficient of determination using broadcasting. Let’s see the implementation

 
 function R²(x, y)
     X = [ones(11) x]
     model = X \ y
     prediction = X * model
     error = y - prediction
     SS_res = sum(v -> v ^ 2, error)
     mean_y = mean(y)
     SS_tot = sum(v -> (v - mean_y) ^ 2, y)
     return 1 - SS_res / SS_tot
 end
  

If we wanted to use broadcasting, we could write:

 
 function R²(x, y)
     X = [ones(11) x]
     model = X \ y
     prediction = X * model
     SS_res = sum((y .- prediction) .^ 2)
     SS_tot = sum((y .- mean(y)) .^ 2)
     return 1 - SS_res / SS_tot
 end
  

As you can see, we have changed formulas for SS_res and SS_tot. In both cases we have used the dot (.) twice. For example, in (y .- prediction) .^ 2 we are broadcasting both subtraction and exponentiation.

By now you know four ways of iteratively applying operations to elements of collections:

• using for loops
• using comprehensions
• using the map function (and other similar higher-order functions that take functions as their arguments)
• using broadcasting

You probably ask yourself in which cases you should use which option. Fortunately, this is mostly a convenience and code readability choice. In your projects use the option that is easiest for you to use, and which results in the most readable code. One of the great features of the Julia language is that all these options are fast. This means that most of the time you do not sacrifice performance by choosing one of them over the other. There are exceptions to this rule in some cases. Apart from the differences that I discussed above, one of the most important exceptions is that if you use the for loop or when you use the map function, then using the Threads module or ThreadsX.jl package you can optionally make the operation take advantage of all cores of your processor. Ability to easily support multi-threaded execution of your code is a feature that distinguishes Julia from R and Python.

That’s all for now. Thanks for reading.