Faster matrix math in R on macOS

If you want faster matrix operations in R on your Mac, you can use Apple's BLAS (Basic Linear Algebra Subprograms) library from their Accelerate framework instead of the library which comes with the R binary that you get from CRAN. (Unless you built R from source yourself.) CRAN recommends against this, saying:

Although fast, it is not under our control and may possibly deliver inaccurate results. [Source: R for Mac OS X FAQ]

So proceed at your own discretion, I suppose? To switch Apple's BLAS:

cd /Library/Frameworks/R.framework/Resources/lib
ln -sf /System/Library/Frameworks/Accelerate.framework/Frameworks/vecLib.framework/Versions/Current/libBLAS.dylib libRblas.dylib

You can verify by running sessionInfo() in R. You should see:

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

P.S. Use ln -sf libRblas.0.dylib libRblas.dylib to revert back to default.

Benchmark 1

library(microbenchmark)
d <- 1e3
x <- matrix(rnorm(d^2), d, d)
microbenchmark(tcrossprod(x), solve(x), svd(x), times = 10L)

With R's default (libRblas.0.dylib):

Unit: milliseconds
          expr       min        lq      mean    median        uq       max
 tcrossprod(x)  276.9307  278.8472  283.9451  282.5854  288.7777  299.6700
      solve(x)  681.5183  696.7888  703.4895  703.6030  712.5728  724.7195
        svd(x) 2672.3574 2692.6614 2701.0442 2698.4279 2722.2180 2730.9727

With Apple's vecLib (libBLAS.dylib):

Unit: milliseconds
          expr        min        lq      mean    median        uq       max
 tcrossprod(x)   8.642364  10.59531  10.60989  10.80781  11.35995  11.62067
      solve(x)  37.205422  41.38051  43.77767  42.34975  42.92999  54.71387
        svd(x) 278.528673 290.26926 304.74482 307.49666 317.56819 325.29980

Benchmark 2

library(microbenchmark)
set.seed(20190604)
n <- 1e3
p <- 1e2
b <- rnorm(p + 1, 0, 10)
x <- matrix(runif(n * p, -10, 10), ncol = p, nrow = n)
y <- cbind(1, x) %*% b + rnorm(n, 0, 2)
microbenchmark(lm(y ~ x))

With R's default (libRblas.0.dylib):

Unit: milliseconds
           expr     min       lq     mean   median       uq      max neval
      lm(y ~ x) 17.0145 20.27288 22.58367 22.16553 24.49879 37.37704   100

With Apple's vecLib (libBLAS.dylib):

Unit: milliseconds
           expr      min       lq    mean   median       uq      max neval
      lm(y ~ x) 5.478486 7.888307 9.210001 8.84202 10.30503 15.80064   100