An AI’s take on happiness, or, the bandsaw problem in machine learning

Computers don’t have feelings. Perhaps they never will – could feelings and emotions really be taught to an artificial intelligence (AI)? Well, maybe. Although many would say that emotions are best described in the abstract language of art and poetry, quantitative research into emotions like happiness is being carried out even as I’m writing this. Earlier this year, a group of researchers from the University of Tokyo, MIT and the Recruit Institute of Technology posted a paper on arXiv describing a database of “100,000 happy moments” – sentences describing moments when people were happy – and how they collected those moments using crowdsourcing.

100,000 happy moments? That’s quite a lot of happiness. Maybe enough to teach an AI what it means to be happy? This is the story of how I tried to do just that, and what a happy AI can tell us about black-box machine learning methods, the importance of knowing when and why things go wrong, and bandsaws. Lots and lots of bandsaws.

 

Teaching a neural network to be happy

So, I had decided to train an artificial intelligence (or more precisely, a neural network) to write messages about happy moments, using the “happy moments” database to learn what such messages can look like (the so-called HappyDB can be downloaded from GitHub). I did this using a recurrent neural network with LSTM units (using this Tensorflow code with some minor modifications). Some years ago, this would have been a major feat, but nowadays it’s a routine task in machine learning (see Botnik or AI Weirdness for lots of great – and funny – examples).

Clearly, an AI with feelings needs to have a name, and the one I came up with was hAIppy (yeah, sorry). Training finished, and it was time for hAIppy to tell me about one of its happy moments. At first, the sentences it wrote were far too long, like this:
I was able to go to the park with my family and I was able to get a new job that I was able to get to see my family and I was able to get a new job that I was able to get to see my family and I was able to get a new job that I was able to get a lot of fun.

After I’d fiddled with the settings to fix that, hAIppy started repeating itself:
I was happy when I was able to spend time with my family and I was able to get a new car for the first time in a long time.
I was happy when I was able to spend time with my family and I was able to get a new car for the first time in a long time.
I was happy when I was able to spend time with my family and I was able to get a new car for the first time in a long time.
I was happy when I was able to spend time with my family and I was able to get a new car for the first time in a long time.

So, I added some rules to keep it from saying the same thing over and over again.

In the end, the results were pretty good. The neural network was able to create novel descriptions of moments of happiness, and most of the time they made perfect sense. Most messages involved family and/or food, so it clearly picked up on the fact that those things are very important for our happiness.

But sometimes the messages made less sense. Sometimes hAIppy decided to become a little philosophical. And sometimes, it went to dark places, as you can see in this example of what a conversation between hAIppy and myself (in which all typewriter font text was generated by hAIppy) could look like:

-When were you happy, hAIppy?
-I was happy that I was able to get out of work I went to the park with my husband and we had a great time.
-That’s nice. And lovely to hear that you’re married! I had no idea. Let’s try this again. When were you happy?
-I watched a movie with my girlfriend and we had a great time.
-Right, right, wonderful… So, tell me about another time you were happy.
-I was happy when I was able to get a new bandsaw, to see the movie and I was very happy when I got to see my grandson.
-A bandsaw…? Well, I’m sure you’ll find some use for that. When else were you happy?
-I was happy when I was able to get a new bandsaw, to create the tape recorder was a very happy moment for me as it is something I was looking forward to achieving and the news made my parents quite happy and proud.
-You built a tape recorder using a bandsaw?! I’m impressed, to say the least..! And good to hear that it made your parents proud. One more time: when were you happy?
-I was happy when I found out that my favorite musician, Juliana Hatfield, is putting me a hug and said: "Daddy, you're back!"
-You’re the father of American musician Juliana Hatfield? Were you making the tape recorder in order to tape her music then? Interesting! I feel like I’m really getting to know you. When else were you happy?
-I was happy when I was able to get a new bandsaw, to see the movie and I was very happy to see the attention directed elsewhere.
-Well, I for one would appreciate it if you directed your attention away from the bandsaw… When else were you happy?
-striven Wars fluff dryer
-What. I don’t even. When else?
-FINISHED IT AND IT LOOKS VERY DIRTY OUTSIDE THE BOX.
-Finished what? Your tape recorder? And why are you screaming at me? Jeez. Please calm down and tell me about when you were happy instead.
-I FEEL VERY HAPPY BECAUSE I LEAVE THE LIFE IN PRESENT BUT HAPPIEST MOMENT WAS WHEN I GAVE BIRTH TO MY SON.
-You leave life in the present? hAIppy, is this a call for help? And why are you still screaming? Can you just tell me about a happy moment, please?
-I do a difficult job at work.
-Join the club, my friend. I’m trying to get a bot to tell me about when it was happy, but it just keeps spewing nonsense at me. Let’s try again.
-The pursuit of happiness is such little task.
-Yeah, trivial, isn’t it…? I’m sure you must’ve been happy loads of times then, so tell me about a time when you were happy.
-We bought a mini zoo.
-Wow. Sounds lovely! And expensive. You must really like animals. Tell me about some other times you were happy.
-I got a free coffee at the store. I was happy when I was a kid. I was happy when I splurged for dessert. I went paddleboarding in an old apartment in NYC. I had sex with my friend whom I haven't seen. My wife came to my fraternity to make breakfast. I played the basketball game because I was at a prestigious festival. I finished a few books at the beginning of the closing of the year. Conducted another theme morning with my mom. We fed this made me happy. At night I got some of my 17 cats at the grocery store I found. I bought a new cat from an online cafe.
-Sounds like you’ve had a busy week! One last time now: tell me about the last time you were happy.
-After I cooked a really good cat.
-You what?! I’m shutting you down, hAIppy.
-I was able to get a new bandsaw, to create my favorite sports team.
-Shutting. You. Down.

There we go. hAIppy likes family, food and cats (and apparently cats=food at times), and has a strange fascination with bandsaws.

 

Why bandsaws matter

All of this was just a bit of fun, of course. Teaching AI about emotions and what’s important to us is in reality going to be a much much greater challenge. But maybe there are a couple of things to learn from hAIppy. While it’s obsession with bandsaws is amusing, it also serves as a reminder that sometimes machine learning models pick up on the wrong things in the data. In fact, I ran three different incarnations of hAIppy using different settings, and all three ended up with a fondness for bandsaws. Not dogs (mentioned in 2,649 moments), not rain (1,496 moments), not TV (814 moments), but bandsaws, mentioned in just one moment out of 100,000.

In this example, it becomes pretty clear to us that something has gone wrong with the network. But in other cases, like when black-box machine learning methods are used to diagnose diseases or determine whether a loan application is approved, problems like this may be difficult to spot. Just like hAIppy greatly overestimates the probability that someone uses the word bandsaw when accounting for why they are happy, a black-box model can overestimate the probability that, say, people with a rare ZIP code default on their loans. The output from the box may look normal on the surface (someone could of course mention bandsaws when explaining why they are happy), but in fact be completely wrong (most people don’t obsess about how great things are now they have a new bandsaw, so neither should our model) – and as with the bandsaw we can only find out that something is wrong when a pattern appears. Mentioning a bandsaw once is perfectly fine, but mentioning it over and over again means that we have a problem.

So when next someone shows you the output from a machine learning model, ask them: where are your bandsaws?

Answering this question is not always easy. Far from it – in fact, in more serious applications, such as medical or financial problems, spotting the bandsaws in the model output is often very difficult. Many machine learning methods, including neural networks, are inherently opaque, and it is hard for us to tell why they made a certain prediction or decision. When those decisions are binary (disease/no disease, loan application approved/rejected), it is virtually impossible to detect when something’s gone wrong just from the output.

There are good reasons to look for bandsaws, too. Sometimes erroneous prediction  are obviously wrong – like when a model predicts that the temperature in Uppsala tomorrow will be 46°C. Others times the errors are much less obvious, and those are the errors we should be the most concerned about. Far too often, we are content with reporting mean squared errors, sensitivity or other metrics of how well a model is does predictions. Those measures are useful, but they don’t tell the whole story. They don’t tell us when we have repeated non-obvious errors (like hAIppy’s bandsaws) in our data, or what those errors are. And therefore, they don’t tell us much about what we should do in order to improve our model.

There are ways to looks for bandsaws, of course. In the above conversation with hAIppy, there are moments where the model suddenly ends up in parts of the sample space (in this case – the set of all words that hAIppy has learned), that it hasn’t explored well enough to make educated guesses of how to continue the sentence. This is what leads to statements like striven Wars fluff dryer. And too often it ends up in a part of the sample space where the best guess for what the next word in the sentence should be is bandsaw. In other applications, this type of problems can be diagnosed by assessing how similar the input to the model is to data in the training set – if a lot of the new data is different from the data that the model has been trained on, there is an increased risk of poor predictions. You can do this using measures like the Mahalanobis distance or by plotting a bubble chart where the first two principal components of your data make up the two axes, bubble colours show whether an observation is from the test or training set and the bubble sizes represent prediction errors. Another interesting option is to run an unsupervised cluster analysis on the test dataset and compare the performance of the model on the different clusters: if you find a cluster for which the performance is significantly worse you may well be on the bandsaws’ trail. And there are many other options for exploring what’s inside of your black box (here, here, here and here).

The point I’m trying to make here is this: you can’t tell from loss functions and common performance measures whether your black-box model has developed an obsession with woodworking machines. Simply looking at those numbers is not enough. To find its bandsaws, you need to run additional experiments on your trained model. Look for repeated mistakes and part of the sample space where it performs poorly. That will tell you where it needs to improve.

 

Transparency is important – for two reasons

Recently, there have been calls (and rightly so!) for accountability and transparency when machine learning is used – if your loan application is rejected by an algorithm, you should have the right to know why.

I agree with the need for this: when machines and algorithms rule our lives, we need insight into how they make their decisions. But for those of use who build the machine learning models there is another reason for trying to open up the black box: it will help us find the bandsaws – the erroneous decisions that look perfectly normal – and thereby help us improve our models. And when you get rid of the bandsaws, you can kick back, relax, and be happy – with family, food… or your very own mini zoo.

Speeding up R computations Pt III: parallelization

In two previous posts, I have written about how you can speed up your R computations either by using strange notation and non-standard functions or by compiling your code. Last year my department bought a 64-core computational server, which allowed me to get some seriously improved performance by using parallelization. This post is about four insights that I’ve gained from my experiments with parallel programming in R. These are, I presume, pretty basic and well-known, but they are things that I didn’t know before and at least the last two of them came as something of a surprise to me.

I’ll start this post by a very brief introduction to parallelization using the doMC package in R, in which I introduce a toy example with matrix inversions. I will however not go into detail about how to use the doMC package.  If you already are familiar with doMC, you’ll probably still want to glance at the matrix inversion functions before you skip ahead to the four points about what I’ve learned from my experiments, as these functions will occur again later in the text.

Parallel programming in R
Modern processors have several cores, which essentially are independent CPUs that can work simultaneously. When your code is parallelized several computations are carried out in parallel on the different cores. This allows you to, for instance, run several steps of a for-loop in parallel. In R it is simple to do this using e.g. the doMC package.

As an example, imagine that you wish to generate and invert several random 200 by 200 matrices and save the inverted matrices in a list. Using a standard for loop, you can do this using the following function:

# A standard function:
myFunc
<-function(B,size)
{
a
<-vector("list", B)
for(i in 1:B)
{
myMatrix
<-matrix(runif(size^2),size,size)
a
[[i]]<-solve(myMatrix)
}
return(a)
}

Using the doMC package, we can write a parallelized version of this function as follows:

library(doMC) # Use doMC package
registerDoMC() # Determine number of cores to use, default = #cores/2

 

# A parallelized function:
myParFunc<-function(B,size)
{
    foreach(i = 1:B) %dopar% 
    {
        myMatrix<-matrix(runif(size^2),size,size)
        solve(myMatrix)
    }
}

 

The parallelized foreach %dopar% loops always returns a list. If you are used to using for-loops in which the output is manually stored in vectors or matrices, this may take some time getting used to. The list structure is actually very useful, and by using unlist, data.frame and Reduce in a clever way you probably won’t have to make a lot of changes in your old code.

Comparing the two functions, we find that the parallelized function is substantially faster:

>  system.time(myFunc(1000,200))
   user  system elapsed 
 19.026   0.148  19.186 
>  system.time(myParFunc(1000,200))
   user  system elapsed 
 24.982   8.448   5.576

The improvement is even greater for larger matrices. Parallelization can save you a lot of time! Unfortunately, there are situations where it actually makes your code slower, as we will see next.

1. Parallelization is only useful when you have computation-heavy loop steps
Each parallelization has an overhead cost in that it takes some time to assign the different loop steps to the different cores. If each step of your loops is quickly executed, the overhead cost might actually be larger than the gain from the parallelization! As an example, if you wish to generate and invert several small matrices, say of size 25 by 25, the overhead cost is too large for parallelization to be of any use:

>  system.time(myFunc(1000,25))
   user  system elapsed 
  0.264   0.004   0.266 
>  system.time(myParFunc(1000,25))
   user  system elapsed 
  0.972   1.480   1.137

In this case, the parallelized function is much slower than the standard unparallelized version. Parallelization won’t magically make all your loops faster!

2. Only parallelize the outer loops
If you are using nested loops, you should probably refuse the urge to replace all your for-loops by parallelized loops. If you have parallelized loops within parallelized loops, you tend to lose much of the gain from the outer parallelization, as overhead costs increase. This is a special case of point 1 above: it is probably the entire inner loop that makes each step of the outer loop computation-heavy, rather than a single step in the inner loop.

Returning to the matrix inversion example, let’s now assume that we wish to use an inner loop to generate a new type of random matrices by multiplying with other matrices with uniformly distributed elements. The unparallelized and parallelized versions of this are as follows:

# No parallelization:
myFunc<-function(B,size)
{
   a<-vector("list", B)
   for(i in 1:B)
   {
    myMatrix<-matrix(runif(size^2),size,size)
    for(j in 1:B)
    {
        myMatrix<-myMatrix %*% matrix(runif(size^2),size,size)
    }
    a[[i]]<-solve(myMatrix)
   }
}

 

# Outer loop parallelized:
myParFunc<-function(B,size)
{
    foreach(i = 1:B) %dopar% 
    {
        myMatrix<-matrix(runif(size^2),size,size)
        for(j in 1:B)
        {
         myMatrix<-myMatrix %*% matrix(runif(size^2),size,size)
        }
        solve(myMatrix)
    }
}

 

# Both loops parallelized:
myParFunc2<-function(B,size)
{
    foreach(i = 1:B) %dopar% 
    {
        myMatrix<-matrix(runif(size^2),size,size)
        myMatrix<-myMatrix %*% Reduce("%*%",foreach(j = 1:B) %dopar%
        {
         matrix(runif(size^2),size,size)
        })
        solve(myMatrix)
    }
}

When the outer loop is parallelized, the performance is much better than when it isn’t:

>  system.time(myFunc(1000,25))
   user  system elapsed 
 89.305   0.012  89.381 
>  system.time(myParFunc(1000,25))
   user  system elapsed 
 46.767   0.780   5.239 

However, parallelizing the inner loop as well slows down the execution:

>   system.time(myParFunc2(1000,25))
   user  system elapsed 
970.809 428.371  51.302 

If you have several levels of nested loops, it may be beneficial to parallelize the two or three outmost loops. How far you should parallelize depends on what is going on inside each loop. The best way to find out how far you should go is to simply try different levels of parallelization.

3. More is not always merrier when it comes to cores
Using a larger number of cores means spending more time on assigning tasks to different cores. When I tested the code for a large-scale simulation study by using different number of cores on our 64-core number cruncher, it turned out that the code executed much faster when I used only 16 cores instead of 32 or 64. Using 16 cores, I might add, was however definitely faster than using 8 cores. Then again, in another simulation, using 32 cores turned out to be faster than using 16. Less is not always merrier either.

An example where a larger number of cores slows down the execution is the following, which is a modified version of myParFunc2 above, in which the multiplication has been replaced by addition to avoid issues with singular matrices:

myParFunc3<-function(B,size)
{
    foreach(i = 1:B) %dopar% 
    {
        myMatrix<-matrix(runif(size^2),size,size)
        myMatrix<-myMatrix %*% Reduce("+",foreach(j = 1:B) %dopar%
        {
         matrix(runif(size^2),size,size)
        })
        solve(myMatrix)
    }
}

 

> registerDoMC(cores=16)
> system.time(myParFunc3(1000,75))
    user   system  elapsed 
2095.571 1340.748  137.397

 

> registerDoMC(cores=32)
> system.time(myParFunc3(1000,75))
    user   system  elapsed 
2767.289 2336.450  104.773

 

> registerDoMC(cores=64)
> system.time(myParFunc3(1000,75))
    user   system  elapsed 
2873.555 3522.984  112.809

In this example, using 32 cores if faster than using either 16 or 64 cores. It pays off to make a small trial in order to investigate how many cores you should use for a particular parallelized function!

4. Make fair comparisons to unparallelized code
The doMC package includes the option of running foreach-loops without parallelization by simply typing %do% instead of %dopar%. This is nice as it lets you compare the performance of parallelized and unparallelized loops without having to rewrite the code too much. It is however important to know that the %do% loop behaves differently than the standard for loop that you probably would use if you weren’t going to parallelize your code. An example is given by the matrix-inversion functions from point 1 above. The %do% function for this is

# A standard (?) function:
myFunc2<-function(B,size)
{
    foreach(i = 1:B) %do% 
    {
        myMatrix<-matrix(runif(size^2),size,size)
        solve(myMatrix)
    }
}

But this function is noticeably slower than the standard function described above:

>  system.time(myFunc(1000,25))
   user  system elapsed 
  0.264   0.004   0.266
>  system.time(myParFunc2(1000,25))
   user  system elapsed 
  0.952   0.044   0.994

When comparing parallelized and unparallelized loops, avoid %doloops, as these make the parallelized loops seem better than they actually are.

I’m honestly not sure why there is such a large discrepancy between the two unparallelized for-loops. I posted a question on this matter on Stack Overflow last year, but to date it has not received any satisfactory answers. If you happen to know wherein the difference lies, I’d much appreciate if you would let me know.

That’s all for now. These are my four insights from playing around with parallelization in R. It is entirely possible that they contradict some well-known principles of parallel programming – I confess that I am almost completely ignorant about the literature on parallelization. Please let me know in the comments if you disagree with any of the points above or if you have counterexamples
!

For more on parallelization in R, check out R bloggers!

doMc is only one of many many packages for parallelization in R.

Speeding up R computations Pt II: compiling

A year ago I wrote a post on speeding up R computations. Some of the tips that I mentioned then have since been made redundant by a single package: compiler. Forget about worrying about curly brackets and whether to write 3*3 or 3^2 – compiler solves those problems for you.

compiler provides a byte code compiler for R, which offers the possibility of compiling your functions before executing them. Among other things, this speeds up those embarrassingly slow for loops that you’ve been using:

> myFunction<-function() { for(i in 1:10000000) { 1*(1+1) } }
> library(compiler)
> myCompiledFunction <- cmpfun(myFunction) # Compiled function

> system.time( myFunction() )
   user  system elapsed 
 10.002   0.014  10.021 
> system.time( myCompiledFunction() )
   user  system elapsed 
  0.692   0.008   0.700 

That’s 14 times faster!

Functions written in C (imported using Rcpp) are still much faster than the compiled byte code, but for those of us who

  • don’t know C,
  • know C but prefer to write code in R,
  • know C but don’t want to rewrite functions that we’ve already written in R,

compiler offers a great way of speeding up computations. It’s included in the recommended R packages since R 2.13 (meaning that it comes with your basic R installation) and since R 2.14 most standard functions are already compiled. If you still are running an older version of R it’s definitely time to update.

Speeding up R computations

The past few days I’ve been going through some R code that I wrote last year, when I was preparing a massive simulation-based power study for some tests for multivariate normality that I’ve been working on. My goal was to reduce the time needed to run the simulation. I wasn’t expecting great improvement, since I’ve always believed that the most common R functions are properly vectorized and optimized for speed. Turns out I was wrong. Very wrong.

The first thing that I did was that I replaced all parentheses ( ) by curly brackets { }. I was inspired to do so by this post (and this, via Xi’Ans Og) over at Radford Neal’s blog. As he pointed out, code that uses parentheses is actually slower than the same code with curly brackets:

> system.time( for(i in 1:1000000) { 1*(1+1) } )
   user  system elapsed 
  1.337   0.005   1.349 
> system.time( for(i in 1:1000000) { 1*{1+1} } )
   user  system elapsed 
  1.072   0.003   1.076 

Similarly, you can compare a*a and a^2:
> system.time( for(i in 1:10000000) 3^2 )
   user  system elapsed 
  5.048   0.028   5.088 
> system.time( for(i in 1:10000000) 3*3 )
   user  system elapsed 
  4.721   0.024   4.748 
So, a^2 is slower than a*a. This made me wonder, are there other built-in R functions that are slower than they ought to be?

One thing that I found very surprising, and frankly rather disturbing, is that mean(x) takes ten times as long to calculate the mean value of the 50 real numbers in the vector x as the “manual” function sum(x)/50:


> x<-rnorm(50)
> system.time(for(i in 1:100000){mean(x)})
   user  system elapsed 
  1.522   0.000   1.523 
> system.time(for(i in 1:100000){sum(x)/length(x)})
   user  system elapsed 
  0.200   0.000   0.200 
> system.time(for(i in 1:100000){sum(x)/50})
   user  system elapsed 
  0.167   0.000   0.167 
> system.time(for(i in 1:100000){ overn<-rep(1/50,50); x%*%overn })
   user  system elapsed 
  0.678   0.000   0.677 
> overn<-rep(1/50,50); system.time(for(i in 1:100000){ x%*%overn })
   user  system elapsed 
  0.164   0.000   0.164 

I guess that the R development core team have been focusing on making R an easy-to-use high level programming language rather than optimizing all functions, but the poor performance of mean is just embarrassing.

Similarly, the var function can be greatly improved upon. Here are some of the many possibilites:


> x <- rnorm(50)
> system.time( for(i in 1:100000) { var(x) } )
   user  system elapsed 
  4.921   0.000   4.925 
> system.time( for(i in 1:100000) { sum((x-mean(x))^2)/{length(x)-1} } )
   user  system elapsed 
  2.322   0.000   2.325 
> system.time( for(i in 1:100000) { {sum(x*x)-sum(x)*sum(x)/length(x)}/{length(x)-1} } )
   user  system elapsed 
  0.736   0.000   0.737 
> system.time( for(i in 1:100000) { {sum(x*x)-sum(x)*sum(x)/50}/49 } )
   user  system elapsed 
  0.618   0.000   0.618 

> system.time( for(i in 1:100000) { sx<-sum(x); {sum(x*x)-sx*sx/50}/49 } )
   user  system elapsed 
  0.567   0.000   0.568

I changed all the uses of mean in my code to “sum/n” instead (and avoided using var entirely) and found that this sped things up quite a bit.
Another trick to speed up your computations is to create the vectors that you wish to change within a loop with the right number of elements. While
a<-NA
for(j in 1:100) a[j]<-j 
works just fine, it is actually quite a bit slower than
a<-rep(NA,100)
for(j in 1:100) a[j]<-j
You could create a in other ways as well of course, for instance by a<-vector(length=100). Here are the numbers:

> system.time( for(i in 1:100000) { a<-NA; for(j in 1:100) a[j]<-j })
   user  system elapsed 
 37.383   0.092  37.482 
> system.time( for(i in 1:100000) { a<-rep(NA,100); for(j in 1:100) a[j]<-j })
   user  system elapsed 
 25.866   0.065  25.936 
> system.time( for(i in 1:100000) { a<-vector(length=100); for(j in 1:100) a[j]<-j })
   user  system elapsed 
 25.517   0.022  25.548

In my case, I’d been a bit sloppy with creating the vectors in my loops in the proper way, so I changed this in my code as well.

In my simulation study, I simulate multivariate random variables, compute some test statistics and use these to estimate the powers of the normality tests against various alternatives. After doing the changes mentioned above, I compared the performance of my old code to that of the new code, for 1000 iterations of the procedure:

> system.time( source(“oldCode.R”) )
   user  system elapsed 
548.045   0.273 548.622 
> system.time( source(“newCode.R”) )

   user  system elapsed 
 93.138   0.002  93.194
The improved code is almost 6 times faster than the old code. When you do ten million or so iterations, that matters. A lot.

In conclusion, it’s definitely possible to speed up your code significantly if you know of the pitfalls of R. I suspect that I’ll be obsessed with finding more pitfalls in the next few weeks, so I’d be thankful for any hints about other weaknesses that R has.

It should probably be mentioned that R is really fast when things are properly vectorized. Last year, a coworker that uses Matlab challenged me to perform a number of matrix computations faster in R than in Matlab. To his great surprise, R won.

As a final remark, I’m now facing a bit of a dilemma. Should I write readable code; a^6; or fast code; a*a*a*a*a*a?

Update: looking to speed up your R computations even more? See my posts on compiling your code and parallelization.