BLAS and LAPACK

Statistics and Linear Algebra

An awful lot of statistics is at its core linear algebra.

For example:

• Linear regession models, find

β̂ = (X^TX)^− 1X^Ty

• Principle component analysis

  • Find T = XW where W is a matrix whose columns are the eigenvectors of X^TX.

  • Often solved via SVD - Let X = UΣW^T then T = UΣ.

Numerical Linear Algebra

Not unique to Statistics, these are the type of problems that come up across all areas of numerical computing.

• Numerical linear algebra  ≠  mathematical linear algebra

• Efficiency and stability of numerical algorithms matter

  • Designing and implementing these algorithms is hard

• Don’t reinvent the wheel - common core linear algebra tools (well defined API)

Modern variants?

Most available BLAS and LAPACK implementations (like R’s defaults) are somewhat dated

• Designed for a single core cpu

• Certain (potentially non-optimal) hard coded defaults (e.g. block size).

Multithreaded alternatives:

• ATLAS - Automatically Tuned Linear Algebra Software

• OpenBLAS - fork of GotoBLAS from TACC at UTexas

• Intel MKL - Math Kernel Library, part of Intel’s commercial compiler tools

• cuBLAS / Magma - hybrid CPU / GPU library from the ICL at UTK

Naming conventions

A BLAS and LAPACK subroutines are named using form pmmaaa where:

• p is a one letter code for the type of data

  • S single precision floating point
  • D double precision floating point
  • C complex single precision floating point
  • Z complex double precision floating point

• mm is a two letter code for the type of matrix expected by the subroutine

• aaa is a one to three letter code denoting the algorithm implemented by subroutine

BLAS Examples

DGEMM

• D - expects data of type double

• GE - expects a general matrix (no useful structure)

• MM - matrix / matrix multiplication operation

STRMV

• S - expects data of type double

• TR - expects a triangular

• MV - matrix / vector multiplication operation

STRMV details

dtrmv(   character   UPLO,
         character   TRANS,
         character   DIAG,
         integer     N,
         double precision, dimension(lda,*)  A,
         integer     LDA,
         double precision, dimension(*)  X,
         integer     INCX 
     )   

STRMV performs one of the matrix-vector operations

x = A x        or         x = A^Tx

where x is an n element vector and A is an n by n unit, or non-unit, upper or lower triangular matrix.

LAPACK Example

DPOTRF

• D - expects data of type double

• PO - positive definite matrix

• TRF - triangular factorization

DPOTRF details

dpotrf(  character   UPLO,
         integer     N,
         double precision, dimension( lda, * )   A,
         integer     LDA,
         integer     INFO 
      )   

DPOTRF computes the Cholesky factorization of a real symmetric positive definite matrix A.

The factorization has the form
A = U^T * U, if UPLO = ’U’, or

A = L * L^T, if UPLO = ’L’,
where U is an upper triangular matrix and L is lower triangular.

Optimization beyond parallelism?

$ /usr/bin/Rscript -e "x=matrix(runif(1000^2),ncol=1000);system.time(x%*%x)"
  user  system elapsed 
 1.258   0.010   1.268

$ /usr/local/bin/Rscript -e "library(OpenBlasThreads);openblas.set.num.threads(1);x=matrix(runif(1000^2),ncol=1000);system.time(x%*%x)"
  user  system elapsed 
 0.123   0.004   0.127 


$ /usr/bin/Rscript -e "x=matrix(runif(5000^2),ncol=5000);system.time(x%*%x)"
   user  system elapsed 
134.687   0.055 134.791 

$ /usr/local/bin/Rscript -e "library(OpenBlasThreads);openblas.set.num.threads(1);x=matrix(runif(5000^2),ncol=5000);system.time(x%*%x)"
   user  system elapsed 
 14.197   0.057  14.250 

OpenBLAS DGEMM Performance

library(OpenBlasThreads)
x=matrix(runif(5000^2),ncol=5000)

sizes = c(100,500,1000,2000,3000,4000,5000)
cores = c(1,2,4,8)

sapply(cores, 
       function(c) 
       {
            openblas.set.num.threads(c)
            sapply(sizes, 
                   function(s) 
                   {
                        y = x[1:s,1:s]
                        system.time(y %*% y)[3]
                   })
       })

OpenBLAS DPOTRF Performance

library(OpenBlasThreads)
library(clusterGeneration)
x=genPositiveDefMat(5000)$Sigma

sizes = c(100,500,1000,2000,3000,4000,5000)
cores = c(1,2,4,8)

sapply(cores, 
       function(c) 
       {
            openblas.set.num.threads(c)
            sapply(sizes, 
                   function(s) 
                   {
                        y = x[1:s,1:s]
                        system.time(chol(y))[3]
                   })
       })

Base BLAS takes ~26 secs for a 5000 x 5000 matrix.