I’m working on optimising some R code written by a researcher at University of Sheffield and its very much a war of attrition! There’s no easily optimisable hotspot and there’s no obvious way to leverage parallelism. Progress is being made by steadily identifying places here and there where we can do a little better. 10% here and 20% there can eventually add up to something worth shouting about.<\/p>\n
One such micro-optimisation we discovered involved multiplying two matrices together where one of them needed to be transposed. Here’s a minimal example.<\/p>\n
#Set random seed for reproducibility\r\nset.seed(3)\r\n\r\n# Generate two random n by n matrices\r\nn = 10\r\na = matrix(runif(n*n,0,1),n,n)\r\nb = matrix(runif(n*n,0,1),n,n)\r\n\r\n# Multiply the matrix a by the transpose of b\r\nc = a %*% t(b)\r\n<\/pre>\nWhen the speed of linear algebra computations are an issue in R, it makes sense to use a version that is linked to a fast implementation of BLAS<\/a> and LAPACK<\/a> and we are already doing that on our HPC system<\/a>.<\/p>\n
Here, I am using version 3.3.3 of Microsoft R Open<\/a> which links to Intel’s MKL <\/a>(an implementation of BLAS and LAPACK) on a Windows laptop.<\/p>\n
In R, there is another way to do the computation c = a %*% t(b)<\/strong>\u00a0 — we can make use of the tcrossprod<\/a> function (There is also a crossprod function for when you want to do t(a) %*% b<\/strong>)<\/p>\n
c_new = tcrossprod(a,b)<\/pre>\nLet’s check for equality<\/p>\n
c_new == c\r\n[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]\r\n[1,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE\r\n[2,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE\r\n[3,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE\r\n[4,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE\r\n[5,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE\r\n[6,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE\r\n[7,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE\r\n[8,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE\r\n[9,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE\r\n[10,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE\r\n<\/pre>\nSometimes, when comparing the two methods you may find that some of those entries are FALSE which may worry you!
\nIf that happens, computing the difference between the two results should convince you that all is OK and that the differences are just because of numerical noise. This happens sometimes when dealing with floating point arithmetic (For example, see https:\/\/www.walkingrandomly.com\/?p=5380<\/a>).<\/p>\nLet’s time the two methods using the microbenchmark package.<\/p>\n
install.packages('microbenchmark')\r\nlibrary(microbenchmark)\r\n<\/pre>\nWe time just the matrix multiplication part of the code above:<\/p>\n
microbenchmark(\r\noriginal = a %*% t(b),\r\ntcrossprod = tcrossprod(a,b)\r\n)\r\n\r\n\r\nUnit: nanoseconds\r\nexpr min lq mean median uq max neval\r\noriginal 2918 3283 3491.312 3283 3647 18599 1000\r\ntcrossprod 365 730 756.278 730 730 10576 1000\r\n<\/pre>\nWe are only saving microseconds here but that’s more than a factor of 4 speed-up in this small matrix case. If that computation is being performed a lot in a tight loop (and for our real application, it was), it can add up to quite a difference.<\/p>\n
As the matrices get bigger, the speed-benefit in percentage terms gets lower but tcrossprod always seems to be the faster method. For example, here are the results for 1000 x 1000 matrices<\/p>\n
#Set random seed for reproducibility\r\nset.seed(3)\r\n\r\n# Generate two random n by n matrices\r\nn = 1000\r\na = matrix(runif(n*n,0,1),n,n)\r\nb = matrix(runif(n*n,0,1),n,n)\r\n\r\nmicrobenchmark(\r\noriginal = a %*% t(b),\r\ntcrossprod = tcrossprod(a,b)\r\n)\r\n\r\nUnit: milliseconds\r\nexpr min lq mean median uq max neval\r\noriginal 18.93015 26.65027 31.55521 29.17599 31.90593 71.95318 100\r\ntcrossprod 13.27372 18.76386 24.12531 21.68015 23.71739 61.65373 100\r\n<\/pre>\n