聊聊 R 和 GPU

注：本文来自寇强的博客，原文请点击此处。

寇强：现为 Indiana University PhD in Informatics。

微博：@没故事的生科男。

这是一直想写几句的一个话题，既然今天有时间就聊一聊。

GPGPU 算是近几年兴起的一个领域，以 CUDA 为代表，在高性能计算方面成果相当多。作为一种相对廉价的高性能解决方案，越来越多的程序员开始加入 GPGPU 阵营。Andrew Ng（就是那个 Machine Learning 公开课的 Andrew）去年在 Google 用造价大约一百万美的集群完成了猫脸识别，而这个月他刚刚宣布他的团队用造价两万美元的 GPU 集群，达到了同样的效果（论文在这儿）。从这里我们也可以大致看出 GPU 在 Machine Learning 方面的潜力。

在这个所谓的 “大数据时代”，大规模机器学习似乎是个必须要解决的问题，而个人观点觉得 Hadoop 这种 MapReduce 平台不适合相对密集型的机器学习。从报道上看，UCBerkely 开发的 Spark 平台在这个方面要远远优于 Hadoop，但 Spark 没玩过，今天还是聊聊以前玩过的 GPU。

GPGPU 的解决方案有不止一个，但由于英伟达集团的大力推广，CUDA 可能是支持最好，也是使用最多的，后面提到的 GPU 也都默认是他家的，所以真正的题目应该是 “R 和 CUDA”。我的测试和开发环境是 ubuntu，后面提到的测试和配置也都是 ubuntu 下面的。

一、GPU 的优势

我们不说 latency 之类的术语，只举个简单的例子。现在全班要去春游，你有一辆保时捷和一辆大巴：保时捷只有四个座位，但半个小时就到了；大巴有 50 个座位，但要一个多小时。为了让全班尽早过去，大巴一定是首选。从计算的角度看，各位的 CPU 就是保时捷，GPU 就是大巴。GPU 每个核心都很弱，但众多的核心还是让 GPU 在并行计算上拥有了相当的优势。另外一点，GPU 有相当的价格优势。单纯从浮点数计算能力来看，300 块左右的 GT430（91.564G）已经接近于 2000 块左右的 i7（107.6G）。

二、GPU 的弱势

简单地讲，不是所有运算都可以并行化，其实这也是并行计算的弱势。但几乎所有矩阵运算都有并行化的可能，所以 Machine Learning 的很多方法移植到 GPU 还是相当有搞头的。

三、用 GPU 的几个 R package

估计多数人都没有精力自己写 CUDA 的代码，我们先来看看和 GPU 有关的几个 package，今天只聊我用过的三个，也就是前三个。

• gputools
• HiPLARM
• magma
• OpenCL
• WideLM
• cudaBayesreg
• permGPU

四、安装和环境配置

现在 CUDA 的安装配置应该已经简单很多了，以 ubuntu 为例，基本只要下载对应的 deb 文件，用 apt-get 就可以完成了，具体请参考 CUDA zone 网站。

安装之后就是配置工作，主要是$PATH（/usr/local/cuda-5.5/bin）、$CUDA_HOME（/usr/local/cuda-5.5）和$LD_LIBRARY_PATH（/usr/local/cuda-5.5/lib64）这几个环境变量的配置。我这里给出的都是默认位置，大家搞清楚自己的 CUDA 安装在哪里就 OK 了。

1. gputools 的安装和示例

gputools 基本上是最具通用性的 package，由 Michigan 大学开发。

这个 package 已经在 CRAN 上了，但直接 install.package 还是可能出错，其实还是环境变量的问题。如果 CUDA_HOME 没有问题的话，检查一下 src 文件夹下的 config.mk，看下面三个变量是不是和自己的一致。如果不一致，应该会报告找不到 R.h 之类的错误。检查过之后就可以 R CMD INSTALL gputools 了。

R_HOME := $(shell R RHOME)
R_INC := $(R_HOME)/include
R_LIB := $(R_HOME)/lib 

这里拿一个矩阵相乘做例子，测试函数如下：

library(gputools)
gpu.matmult <- function(n) {
    A <- matrix(runif(n * n), n ,n)
    B <- matrix(runif(n * n), n ,n)
    tic <- Sys.time()
    C <- A %*% B
    toc <- Sys.time()
    comp.time <- toc - tic
    cat("CPU: ", comp.time, "\n")
    tic <- Sys.time()
    C <- gpuMatMult(A, B)
    toc <- Sys.time()
    comp.time <- toc - tic
    cat("GPU: ", comp.time, "\n")
}

可以明显地看到，在维度比较低的时候，GPU 没有任何优势。

gpu.matmult(5)

## CPU:  0.0001199245 
## GPU:  0.2105169 

gpu.matmult(50)

## CPU:  0.000446558 
## GPU:  0.003529072 

gpu.matmult(500)

## CPU:  0.07863498 
## GPU:  0.02003336 

开始有优势了！

gpu.matmult(1000)

## CPU:  0.7417495 
## GPU:  0.09884238 

gpu.matmult(2000)

## CPU:  5.727753 
## GPU:  0.7211812 

2. rpud 的安装和示例

rpud 是著名的 R Tutorial 网站开发的，下载界面提供了三个 package，RPUD、RPUDPLUS 和 RPUSVM。其中只有 rpud 是开源的，后两个只提供了编译好的文件。

rpud 的安装就容易多了，添加一个 R_LIBS_USER 环境变量，之后直接安装即可。

拿 Matrix Distance 举个例子：

test.data <- function(dim, num, seed = 17) {
    set.seed(seed)
    matrix(rnorm(dim * num), nrow = num)
}
m <- test.data(120, 4500)
system.time(dist(m))

##    user  system elapsed 
##  13.944   0.016  13.977

library(rpud)

## Rpud 0.3.4
## http://www.r-tutor.com
## Copyright (C) 2010-2013 Chi Yau. All Rights Reserved.
## Rpud is licensed under GNU GPL v3. There is absolutely NO warranty.

system.time(rpuDist(m))

##    user  system elapsed 
##   0.452   0.248   0.702

加速效果还不错，但不开源就着实让人不爽了。

3. HiPLAR 的安装和示例

HiPLAR 是 High Performance Linear Algebra in R 的缩写，这个 package 的配置略复杂，因为调用的 library 略多。不过好在提供了 installer，一般也不会出错。这里直接拿官方例子说话吧。

library(Matrix);
n <- 8192;
X <- Hilbert(n);
A <- nearPD(X);
system.time(B <- chol(A$mat));

## user  system elapsed  
## 97.990   0.356  98.591

library(HiPLARM)
system.time(B <- chol(A$mat));

## user  system elapsed 
## 1.012   0.316   1.337

五、调用 CUDA

用了各种 package，但很多时候还是要自己写 CUDA 代码，这个也是我正在做的。

这里用的是 Matloff 老爷子的 Programming on Parallel Machines 里的例子。具体的 CUDA 语法这里就不聊了，我入门是 Udacity 上的视频，这里也推荐给大家。

CUDA 是基于 C 的，所以调用 CUDA 和调用 C 没有多少不同。调用 C 的时候，是把 C 文件用 R CMD SHLIB 编译成 so，之后用 R 调用。稍微注意一下 R CMD SHLIB 的输出，其实就什么都明白了。

$ R CMD SHLIB sd.c
gcc -std=gnu99 -I/usr/share/R/include -DNDEBUG -fpic  -O3 -pipe  -g  -c sd.c -o sd.o
gcc -std=gnu99 -shared -o sd.so sd.o -L/usr/lib/R/lib -lR

R CMD SHLIB 只是自动调用了下面那两行，对于 CUDA，我们手动写就行了。

CUDA 代码如下：

#include<cuda.h>;
#include<stdio.h>;

extern "C" void meanout(int *hm, int *nrc, double *meanmut);

__device__ void findpair(int tn, int n, int *pair)
{
    int sum = 0, oldsum = 0, i;
    for(i = 0; ; i++){
        sum += n - i - 1;
        if(tn <= sum - 1){
	pair[0] = i;
	pair[1] = tn - oldsum + i + 1;
	return;
	}

        oldsum = sum;
    }

}

__global__ void proc1pair(int *m, int *tot, int n)
{
    int pair[2];
    findpair(threadIdx.x, n, pair);
    int sum = 0;
    int startrowa = pair[0], startrowb = pair[1];
    for (int k = 0 ; k < n ; k++)
        sum += m[startrowa + n*k]*m[startrowb + n*k];
    atomicAdd(tot, sum);
}

void meanout(int *hm, int *nrc, double *meanmut)
{
    int n = *nrc, msize = n*n*sizeof(int);
    int *dm, htot, *dtot;
    cudaMalloc((void **)&dm, msize);
    cudaMemcpy(dm, hm, msize, cudaMemcpyHostToDevice);
    htot = 0;
    cudaMalloc((void **)&dtot, sizeof(int));
    cudaMemcpy(dtot, &htot, sizeof(int), cudaMemcpyHostToDevice);

    dim3 dimGrid(1, 1);
    int npairs = n*(n - 1)/2;
    dim3 dimBlock(npairs, 1, 1);
    proc1pair<<<dimGrid, dimBlock>>>(dm, dtot, n);
    cudaThreadSynchronize();
    cudaMemcpy(&htot, dtot, sizeof(int), cudaMemcpyDeviceToHost);
    *meanmut = htot/double(npairs);
    cudaFree(dm);
    cudaFree(dtot);
}

编译选项如下：

$ nvcc -g -G -I/usr/local/cuda/include -Xcompiler "-I/usr/share/R/include -fpic" -c mutlinksforr.cu -o mutlink.o -arch=sm_11
$ nvcc -shared -Xlinker "-L/usr/lib/R/lib -lR" -L/usr/local/cuda/lib mutlink.o -o meanlinks.so

R 里的调用和输出

dyn.load("meanlinks.so")
m <- rbind(c(0, 1, 1, 1), c(1, 0, 0, 1), c(1, 0, 0, 1), c(1, 1, 1, 0))
ma <- rbind(c(0, 1, 0), c(1, 0, 0), c(1, 0, 0))
.C("meanout", as.integer(m), as.integer(4), mo = double(1))

## [[1]]
##  [1] 0 1 1 1 1 0 0 1 1 0 0 1 1 1 1 0
## 
## [[2]]
## [1] 4
## 
## $mo
## [1] 1.333

今晚有时间，所以就大致聊了聊 R 和 GPU，准确地说是 R 和 CUDA。利用 GPU 做 machine learning 是我现在除去实验室项目之外，最大的兴趣所在。现在先从 Back-Propagation 开始，也许明年会有一个基于 GPU 的 machine learning 的 package 出来，但不保证开发进度呀。（如果真的发布 0.1 版本，不知道中国 R 语言大会能不能混个演讲。）

特别提一点，这是我笔记本上的测试结果，显卡是 GT 720M，这个并不是专业的计算卡，如果动用 Tesla，效果应该明显得多。刚刚拿到学校的 GPU 集群的帐号，哪天可以试一试。

寇强 http://weibo.com/thirdwing?is_all=1，Rcpp 核心团队成员。本科就读于中山大学，现为印第安纳大学博士在读。

敬告各位友媒，如需转载，请与统计之都小编联系（直接留言或发至邮箱：[email protected]），获准转载的请在显著位置注明作者和出处（转载自：统计之都），并在文章结尾处附上统计之都微信二维码。

← COS 每周精选: 21 世纪的贝叶斯定理？ 第六届中国 R 语言会议（上海会场）通知 →

发表 / 查看评论