I am currently developing my first R package which confronted me a lot with the question: “How can I speed up my code?”.

I did some “research” and read a lot of general articles about speeding up code, but also a few posts specifically about speeding up R code. While I mostly got the main points, I always found the example use cases slightly contrived. So I decided to write a little something about one of my use cases which includes many points that I think are important when trying to speed up your code. The method we want to implement in this post is the so called Bradley Terry Model. If you do not care about the theoretical part, you can jump directly to the implementation section.

### The Bradley Terry Model

The Bradley Terry model is a fairly popular probabilistic model for paired comparisons in different contexts. In sports, for example, paired comparisons of course mean matches between teams or individual players. Given a pair of teams/players i and j, we can estimate a probability that “i beats j”, or i>j, as

$$P(i>j)=\frac{p_i}{p_i+p_j}.$$

The values and represent a skill rating of teams or players that were somehow (we come to that later) estimated before. If we know the skills, thus represents the probability that i will win a match against j. Conversely,

$$P(j>i)=\frac{p_j}{p_i+p_j}=1-P(i>j).$$

This of course only holds if there are no ties involved, as for example in tennis, which will be our use case for the method.

#### Estimating skills

The player or team skills can be estimated from past games with the log-likelihood function

$$l(p)=\sum\limits_{i=1}^n\sum\limits_{i=1}^n \left(w_{ij}\log(p_i)-w_{ij}\log(p_i+p_j)\right),$$

where is the the number of times i has beaten j. The log-likelihood can be maximized iteratively with

$$p_i^{(k+1)}=W_i\left(\sum\limits_{j \neq i} \frac{N_{ij}}{p_i^{(k)}+p_j^{(k)}}\right)^{-1},$$

where is the number of wins by i and is the number of games between i and j. After each iteration, the skill values are scaled to sum to one. The procedure is stopped as soon as

$$\vert p^{(k+1)}-p^{(k)}\vert \leq \epsilon$$

for some .

The alert mathematician might notice that there is a slight problem with the iterative procedure. If there is a team/player i that does not win any game, the procedure fails. A safety net is to give every team a teeny tiny victory against every other team. In math terms, we make irreducible such that the iterative procedure converges.

Once the process converges, we have an estimate of player/team skills that can be used to form a ranking.

### A first implementation

The iterative process for is easily implemented in R.

BTm.1=function(matches,max.iter=100,tol=10^-8){ #matches should have two columns, where the first column should always be the winner g=graph_from_edgelist(as.matrix(matches),directed=T) E(g)$weight=1 g=igraph::simplify(g,edge.attr.comb = "sum") A=get.adjacency(g,attr = "weight",sparse=F) n=dim(A)[1] # sparse correction if(!is.connected(g,mode = "strong")){ eps=1/n A=A+matrix(eps,n,n)-diag(eps,n) } #no games between players N=A+t(A) #no of wins for each player W=rowSums(A) #initialisation w_0=rep(1/n,n) w_old=rep(n,n) iter=0 while(iter<=max.iter & sum((w_old-w_0)^2)>tol){ iter=iter+1 w_old=w_0 for(i in 1:n){ tmp=0 for(j in 1:n){ if(i!=j){ tmp=tmp+N[i,j]/(w_0[i]+w_0[j]) } } w_0[i]=W[i]*tmp^(-1) } w_0=w_0/sum(w_0) } return(data.frame(player=V(g)$name,skill=w_0)) }

The construction of the game matrix `N`

might seem a bit awkward but I am a network analyst, so constructing the game graph feels natural and was the quickest I could come up with. After constructing the matrix we check if it is irreducible. If not we add a small `eps`

to every entry. The heart of the function is the while loop where the iterative process is carried out.

To benchmark the function, I will use some tennis data from a previous post gathered from here and here. For the beginning we use all ATP matches played in 2016.

atp=read_csv("atp.csv") data1<-atp %>% dplyr::filter(date>as.Date("2016-01-01")) %>% select(winner,loser) glimpse(data1) ## Observations: 2,626 ## Variables: 2 ## $ winner <chr> "Grigor Dimitrov", "Denis Kudla", "Tobias Kamke", "Hyeo... ## $ loser <chr> "Gilles Simon", "John Patrick Smith", "Benjamin Mitchel...

How fast is our implementation?

system.time(BTm.1(data1)) ## user system elapsed ## 11.452 0.000 11.371

11 seconds seems a bit slow for only 2626 observations. Let’s see what we can do about that.

### Maybe there is a Package for it?

Actually this is a question you should ask yourself before implementing a potentially complex function in R. Chances are somewhere, someone has already done that and probably better than you (no offense!). For our problem, there actually exists the `BradleyTerry2`

package.

library(BradleyTerry2) library(microbenchmark) #prepare input data for the Package function data2<-data.frame(player1=factor(data1$winner,levels=unique(c(data1$winner,data1$loser))), player2=factor(data1$loser,levels=unique(c(data1$winner,data1$loser))), outcome=rep(1,nrow(data1))) res<-microbenchmark(BTm.1=BTm.1(data1), package=BTm(data2$outcome,data2$player1,data2$player2), times=10) print(res) ## Unit: seconds ## expr min lq mean median uq max neval ## BTm.1 10.844524 10.93058 11.274145 11.181254 11.40294 12.486450 10 ## package 4.604094 4.66283 4.823965 4.753485 4.99618 5.274492 10 res %>% group_by(expr) %>% dplyr::summarise(mean.t=mean(time)/10^9) %>% ggplot(aes(x=expr,y=mean.t))+ geom_col()+ labs(x="Expression",y="Mean Time in seconds")

The code of the package is almost 3x faster, which was to be expected.

Using functions from existing packages is usually a good idea. Why reinventing the wheel?

But sometimes you might reconsider this, especially when you want to include a function into your own R package.

I personally do not like my package to depend on to many others since I don’t want to force users to install a bulk of packages, that they might not need otherwise. Also, other packages might be using different data structures forcing you to use them as well and you will find yourself transforming data back and forth, slowing down your code again.

If you decide to use your own implementation, you can use existing packages as benchmark for speed and also for accuracy.

res.1=BTm.1(data1) res.2=BTm(data2$outcome,data2$player1,data2$player2) #top five ATP players of 2016 according to BTm.1 res.1$player[order(res.1$skill,decreasing=T)[1:5]] ## [1] "Andy Murray" "Novak Djokovic" "Milos Raonic" "Roger Federer" ## [5] "Kei Nishikori" #top five ATP players of 2016 according to BTm names(sort(res.2$coefficients,decreasing=T))[1:5] ## [1] "..Andy Murray" "..Novak Djokovic" "..Milos Raonic" ## [4] "..Roger Federer" "..Kei Nishikori"

While this is certainly not a very robust test to check if the two functions produce the same outcome, it shows that they at least agree upon the 5 best players in 2016. I am not big tennis connoisseur, but I think this Top 5 looks very reasonable. If you plan on using a self-implemented function more frequently, you should definitely make sure that it is correct.

### Can we speed up our implementation?

Now that we decided to use our own implementation, it is time to think about speeding up the code.

A very easy way to do so code is using a byte code compiler. Explaining the details of this is out of scope for this blog post. A post about this topic can be found here.

BTm.1.jit=compiler::cmpfun(BTm.1) res<-microbenchmark(BTm.1=BTm.1(data1), BTm.1.jit=BTm.1.jit(data1), package=BTm(data2$outcome,data2$player1,data2$player2), times=10) print(res) ## Unit: seconds ## expr min lq mean median uq max ## BTm.1 10.809827 10.908028 10.925320 10.928710 10.958387 10.981341 ## BTm.1.jit 4.468739 4.489647 4.553959 4.538044 4.605068 4.697536 ## package 4.512514 4.591104 4.607991 4.615787 4.634948 4.725338 ## neval ## 10 ## 10 ## 10 res %>% group_by(expr) %>% dplyr::summarise(mean.t=median(time)/10^9) %>% ggplot(aes(x=expr,y=mean.t))+ geom_col()+ labs(x="Expression",y="Median Time in seconds")

With just one line of code we could speed up our function so much, that we now match the runtime of the dedicated package!

### Can we speed things up even more?

Although using for loops is easy and readable, it mostly has a tremendous effect on runtime in R.

This stackoverflow thread discusses some loops in R in some detail.

The key to fast R code is vectorization. The downside of vectorization is that you have to invest some time into it and it requires some knowledge in matrix algebra. Carefully investigating our code we might find that the line `tmp=tmp+N[i,j]/(w_0[i]+w_0[j])`

looks suspiciously like something we could rewrite in terms of matrices. The denominator can be expressed as a matrix `W_0`

where `W_0[i,j]=w_0[i]+w_0[j]`

. Calculating this matrix can either be done as before in a for loop, yet a more efficient way is using the `outer`

function.

Here is a small example.

#for loops f.for=function(a){ n=length(a) A=matrix(0,n,n) for(i in 1:n){ for(j in 1:n){ A[i,j]=a[i]+a[j] } } return(A) } #using outer f.outer=function(a) {return(outer(a,a,"+"))} #are the results identical? n=25 a=runif(n) identical(f.for(a),f.outer(a)) ## [1] TRUE #runtime n=2500 a=runif(n) res<-microbenchmark::microbenchmark(f.for=f.for(a),f.outer=f.outer(a),times=10) print(res) ## Unit: milliseconds ## expr min lq mean median uq max ## f.for 8314.93612 8440.52418 8480.5129 8462.018 8571.8711 8616.8158 ## f.outer 32.84932 33.16897 128.5399 167.679 172.0675 300.8795 ## neval ## 10 ## 10 res %>% group_by(expr) %>% dplyr::summarise(mean.t=median(time)/10^9) %>% ggplot(aes(x=expr,y=mean.t))+ geom_col()+ labs(x="Expression",y="Median Time in seconds")

The difference between the for loops and the outer is enormous. Not only does it cut the required lines of code to a bare minimum, but it is also 60x faster in this case!

Using the outer function in our implementation allows us to reduce each iteration to a simply `rowSums`

call.

BTm.2=function(matches,max.iter=100,tol=10^-8){ #matches should have two columns, where the first column should always be the winner sparse.corrected=F g=graph_from_edgelist(as.matrix(matches),directed=T) E(g)$weight=1 g=igraph::simplify(g,edge.attr.comb = "sum") A=get.adjacency(g,attr = "weight",sparse=F) n=dim(A)[1] # sparse correction if(!is.connected(g,mode = "strong")){ eps=1/n A=A+matrix(eps,n,n)-diag(eps,n) sparse.corrected=T } #no games between players N=A+t(A) #no of wins for each player W=rowSums(A) #initialisation w_0=rep(1/n,n) w_old=rep(n,n) iter=0 while(iter<=max.iter & sqrt(sum((w_old-w_0)^2))>tol){ iter=iter+1 w_old=w_0 w_0=W*rowSums(N/outer(w_0,w_0,"+"))^(-1) w_0=w_0/sum(w_0) } return(data.frame(player=V(g)$name,skill=w_0)) }

Notice how the vectorized code reduced the size of the while loop significantly. Admittedly, it also made it a bit harder to understand. I have learned (the hard way) that when you vectorize your code, always keep the non-vectorized code as a reference.

How much runtime do we gain?

res<-microbenchmark(BTm.1.jit=BTm.1.jit(data1), BTm.2=BTm.2(data1), package=BTm(data2$outcome,data2$player1,data2$player2), times=10) print(res) ## Unit: milliseconds ## expr min lq mean median uq max ## BTm.1.jit 4454.3048 4473.1272 4515.7012 4485.3681 4595.2181 4615.3122 ## BTm.2 431.4755 445.0542 480.6702 458.8493 477.7227 595.9364 ## package 4305.5614 4517.1449 4594.6055 4609.5872 4632.7367 4863.8251 ## neval ## 10 ## 10 ## 10 res %>% group_by(expr) %>% dplyr::summarise(mean.t=median(time)/10^9) %>% ggplot(aes(x=expr,y=mean.t))+ geom_col()+ labs(x="Expression",y="Median Time in seconds")

The speed up of the vectorized code is extraordinary. It is even several magnitudes faster than the code from the `BradleyTerry2`

package. Most likely due to the fact that the `BTm`

function is more generic and also produces a variety of diagnostics that we completely ignored. The code for `BTM.2`

is highly specific but exactly what we wanted.

**Disclaimer**: It might well be that the `Btm`

function could be speed up when using proper options, but I did not spent to much time looking into the details of that function.

### Scaling it up

A speed up is always nice, especially when it is as large as in our case. But a speed up for a single run is not all we need. It is also important how the runtime scales with the size of the input.

Below I increased the time window for the considered matches used in `BTm.2`

gradually from 2016 back to 2000, gradually including more matches.

years=2016:2000 res.scale=data.frame(runtime=rep(0,length(years)), observations=rep(0,length(years))) k=0 for(y in years){ data3<-atp %>% dplyr::filter(lubridate::year(date)>=y) %>% select(winner,loser) k=k+1 tmp<-microbenchmark(BTm.2=BTm.2(data3),times=5) res.scale$runtime[k]=median(tmp$time)/10^9 res.scale$observations[k]=nrow(data3) } ggplot(res.scale,aes(x=observations,y=runtime))+ geom_point()+ geom_line()+ labs(x="Observations",y="Median time in seconds")

Not sure what is going on in the middle (I should have maybe done more than five runs) but apart from that I think the code scales quite reasonably.

### Wrap up

The key to speeding up R code without a doubt is to vectorize your code. This will hardly be any news for more experienced R users. Although it was also nothing new for me, things get a bit more real when you try to develop a package. It might be no problem if you wait for results a couple of seconds (or minutes, if you are patient), but users might have bigger inputs then you have, so scalability becomes an issue too.

I have always been a huge fan of the `outer`

function since it appeals to my inner mathematician. I will dedicate the next post to it in more detail, since you can do impressive things with it.

I did, on purpose, not tackle the option of implementing functions in C++ which will also be left for the next post. But be sure that even vectorized R code can be magnitudes slower than a naive C/C++ implementation.

The code of this post can be found on Github.

By the way, for those who are interested in the ranking of tennis players from 2000-2016. Here is what we get.

Rank | Player | Skill |
---|---|---|

1 | Novak Djokovic | 0.0154895 |

2 | Rafael Nadal | 0.0131991 |

3 | Roger Federer | 0.0128205 |

4 | Andy Murray | 0.0097758 |

5 | Andre Agassi | 0.0061090 |

6 | Andy Roddick | 0.0060375 |

7 | Juan Martin Del Potro | 0.0059108 |

8 | Jo Wilfried Tsonga | 0.0051155 |

9 | Lleyton Hewitt | 0.0047013 |

10 | Milos Raonic | 0.0044824 |