对统计算法进行分析所带来工作负载可能需要在单独一台计算机上运行数小时甚至数天时间。为了解决这个难题，众多统计学家与数据科学家利用R统计计算环境编写出了复杂的模拟与统计分析方案。然而这些成果通常也需要经过漫长的运行处理时间。为了节约R语言程序员们耗费在等待运算结果上的宝贵时间，充分利用计算机制的并行性以及硬件性能自然成了一项极具现实意义的工作。
在之前的一篇文章中，我们曾经探讨过以并行方式处理长时间运行的R模型的价值，同时演示了如何利用多计算核心与多节点并行机制削减运行时长。在今天的文章中，我们将继续探索R语言的另外一种并行计算途径，即利用通用图形处理单元（简称GPU）当中的处理核心显著加快R环境下的聚类算法处理速度。目前普及程度最高的GPU计算方案要数英伟达公司推出的Tesla系列产品了。单一Tesla K40 GPU拥有2880个集成核心、12 GB显存以及高达每秒288 GB的传输带宽，这使其能够在一秒钟之内完成5万亿次浮点运算。
本文所列出的示例由Chi Yau发布于r-tutor.com网站的成果。Chi是CRAN开源rpud软件包以及rpudplus的作者，这些R库旨在帮助开发人员更轻松地发挥GPU的强大处理能力，同时又无需直接在CUDA C++环境中编写代码。要深入了解R语言以及利用GPU实现并行编程，大家可以点击此处下载Chi所撰写的电子书。为了便于说明，今天的文章将把重点放在与距离计算及分级聚类方面，不过大家也可以利用rpud软件包对多种应用程序进行加速。
在R语言环境下进行分级聚类
聚类分析，或者简称聚类，是指将多个彼此相似（基于特定指标）的对象汇聚在同一分组当中并加以分析的过程。聚类分析是并发编程领域的一类常见问题。在之前的文章当中，我们曾经提到过利用k-means实现聚类分析。不过这一次，我们将着眼于在R语言环境下利用hclust实现聚类分析——hclust这一函数能够帮助用户轻松根据观察得出的差异创建出系统树图（即图一所示的树状图）。这类分析工作适用于多种应用程序，包括癌症研究以及针对金融数据的时间序列分析等。
图一：利用R环境下聚类分析创建得出的系统树图。
与此前演示过的k-means示例相似，以分级方式对观察结果进行分组取决于我们是否能够对不同观察结果间的差异（或者距离）进行量化。也就是说对不同观察结果对间的欧氏几何距离加以计算（大家可以将其视为被推广至更多维度当中的勾股定理）。ChiYau曾在他的两篇博文——《GPU实现之距离数组》以及《聚类分析》——当中对此进行过阐述，因此我们就不在这里陈述具体细节了。
R语言的hclust函数能够接受由此前计算得出的观察结果间距离所构成的数组。而R语言中的dist函数则负责计算数据集内不同行间的差异，且支持包括欧氏几何距离（默认情况）在内的多种方法。如果我们现有一组M观察值（行），而每一行都拥有N属性（列），那么要获取每一项距离计算结果，我们需要计算在N维空间当中不同观察值之间的向量长度。所有行之间的离散距离计算结果共有
个。因此，实际计算规模与观察值数量的平方呈正比：对于10个观察值，我们需要进行45次距离计算，而对1000个观察值，我们则需要进行4950次距离计算。而在面对10万个观察值时，我们需要完成的距离计算数量则达到4999950000（接近50亿）次。如大家所见，对于大型数据集来说，距离计算往往会带来高昂的计算成本。
部署GPU环境以完成测试
在我们开始运行应用程序之前，首先需要访问一套拥有GPU的系统。幸运的是，现在我们能够以合理的价格租用到一台拥有GPU的设备。在筹备本文的过程中，我用了两个小时通过租来的设备尝试TeraprocR分析集群即服务。这项服务采用Amazon EC2实例，且租赁的整体成本为每小时1.3美元——这可比购买并安装自有设备便宜多了！我之所以只用了这么点时间就完成了体验过程，是因为该集群的整个安装过程完全由Teraproc服务自动完成。Teraproc的R集群即服务提供CUDA、R、R Studio以及其它需要预先安装且随时能够使用的软件组件。OpenLava与NFS同样以自动化方式配置完成，这让我能够将集群延伸至多个具备GPU的计算节点当中，并可以选择使用Amazon精确实例方案来降低使用成本。
如后文中的图片所示，我利用Amazon g2.2xlarge设备类型在Teraproc.com上部署了一套单节点集群。我本来可以通过AmazonEC2控制台自行安装g2.2xlarge实例的。但这样一来，我还需要自行安装R、R Studio并对整套环境进行配置，这无疑会带来更加高昂的时间与经济成本。大家可以点击此处查看Teraproc R分析集群即服务网站，从中学习如何在不同节点类型（包括免费设备）上亲手建立R集群。如果大家已经拥有Amazon EC2账户，则可以在不到五分钟的时间内完成集群的设置工作。
我使用的这套g2.2xlarge机器实例采用基于SandyBridge架构的八核心/vCPU至强E5-2670处理器、内存容量约为15 GB、配备固态硬盘以及英伟达GRID K520GPU。这套按需实例的每小时使用成本为0.65美元。其中英伟达GRIDK520拥有两块GK104图形处理器，每一块各自拥有1536个计算核心以及8 GB显存。
运行示例
首先，我们利用Teraproc.com的R即集群服务对R环境进行配置，从而确保我们选择了正确的机器类型（G22xlarge）并安装了一套单节点集群，如图一所示。该服务配合R Studio自动部署了一套单节点集群，并提供了一条URL用于访问该R Studio环境。
图二：迈出利用TeraprocR分析集群即服务构建GPU加速型R Studio集群的第一步。
利用RStudio当中的shell功能（位于Tools菜单当中），我可以运行一条操作系统命令来确保该GPU切实存在于该机器当中。
gordsissons@ip-10-0-93-199:~$ lspci | grep -i nvidia
00:03.0 VGA compatible controller: NVIDIA Corporation GK104GL [GRID K520] (rev a1)
要使用rpud软件包访问各类GPU功能，我们需要首先安装该软件包。在R Studio当中运行以下命令，即可将来自CRAN的rpud安装到当前环境中。
> install.packages(“rpud”)
接下来，使用library命令来访问各rpud功能。
> library(rpud)
Rpudplus 0.5.0
http://www.r-tutor.com
Copyright (C) 2010-2015 Chi Yau. All Rights Reserved.
Rpudplus is free for academic use only. There is absolutely NO warranty.
如果一切正常，那么我们应该能够通过在R命令提示符中调用rpuGetDevice来查看到Amazon实例内的GPU。
> rpuGetDevice()
GRID K520 GPU
[1] 0
下面列出的是一款示例R语言程序，其作用是比较使用GPU加速机制与不使用两种情况下聚类分析算法的性能水平。我们要做的第一步是创建一个合适的数据集，其中包含我们能够控制的观察值（即行）数量以及每个观察值的维度数量（也就是列）。其中test.data函数会根据行与列的数量返回一个随机值数组。
run_cpu函数负责利用R的dist函数计算所有不同观察值（也就是行）间的距离，而后针对计算得出并被保存在dcpu内的距离值运行R语言原生hclust函数，从而创建出一套系统树图。而run_gpu函数执行的是同样的计算流程，只不过使用的是rpud软件包当中针对GPU进行过优化的dist与hclust版本（即rpuDist与rpuHclust<em>）。
这套R脚本会调用test.data以创建出一个特定大小的数组m，而后衡量并显示同时利用CPU与GPU函数来创建一个系统聚类时所需要的时间。
library("rpud")#
# function to populate a data matrix
#
test.data <- function(dim, num, seed=10) {
set.seed(seed)
matrix(rnorm(dim * num), nrow=num)
}
run_cpu <- function(matrix) {
dcpu <- dist(matrix)
hclust(dcpu)
}
run_gpu <- function(matrix) {
dgpu <- rpuDist(matrix)
rpuHclust(dgpu)
}
#
# create a matrix with 20,000 observations each with 100 data elements
#
m <- test.data(100, 20000)
#
# Run dist and hclust to calculate hierarchical clusters using CPU
#
print("Calculating hclust with Sandy Bridge CPU")
print(system.time(cpuhclust <-run_cpu(m)))
#
# Run dist and hclust to calculate hierarchical clusters using GPU
#
print("Calculating hclust with NVIDIA K520 GPU")
print(system.time(gpuhclust <- run_gpu(m)))
运行该脚本会得出以下结果：
>source('~/examples/rgpu_hclust.R')[1] "Calculating hclust with Sandy Bridge CPU"
user system elapsed
294.760 0.746 295.314
[1] "Calculating hclust with NVIDIA K520 GPU"
user system elapsed
19.285 3.160 22.431
为了探索GPU与CPU之间的加速能力差异，我们针对数据集内的不同行数运行该脚本并对结果进行整理。距离计算任务在GPU上以高度并发方式进行，而大部分大部分GPU优化型hclust的计算任务则运行在CPU之上。有鉴于此，当数据集规模持续增大且距离计算所占比例由此提高时，GPU加速效果也变得愈发突出。
行数 | 维数 | 总体元素数量 | # 距离计算 | CPU时间（秒） | GPU时间（秒） | 加速效果 |
1,000 | 100 | 100,000 | 1,998,000 | 0.50 | 0.04 | 11.8 |
2,000 | 100 | 200,000 | 7,996,000 | 2.06 | 0.17 | 12.1 |
5,000 | 100 | 500,000 | 49,990,000 | 13.42 | 1.17 | 11.5 |
10,000 | 100 | 1,000,000 | 199,980,000 | 59.83 | 5.03 | 11.9 |
15,000 | 100 | 1,500,000 | 449,970,000 | 141.15 | 11.61 | 12.2 |
20,000 | 100 | 2,000,000 | 799,960,000 | 295.31 | 22.43 | 13.2 |
对运行时间进行逐一比较，我们会看到运行多步骤计算时GPU的速度表现相较单独CPU可提升十倍以上。
我们对该系统聚类的分析结果可通过R语言的plot命令显示为一套系统树图，如图一所示。
> plot(gpuhclust,hang = -1)
总结
我们的研究结果清楚地表明，在GPU上运行此类分析可谓意义重大。我们不仅能够以仅为原本十分之一的时耗完成同样的计算任务，更重要的是可以借此大幅降低完成工作所需要的资源成本。我们可以利用这些时间与资源进行更为深入的分析或者探索更多场景。在Teraproc服务的帮助下，我们能够将GPU的强大计算能力更轻松地交付到R程序员手中，帮助这些一直被困在原地的朋友顺利投向配备GPU节点的怀抱。
