HPC 与 R：基础知识
我们将从一点历史和一些基础知识开始，这些知识是任何接触 R 并关心性能的人都应该了解的。
点击此处注册并下载免费的 Intel® Parallel Studio XE 30天试用版
Drew Schmidt，研究生研究助理，田纳西大学诺克斯维尔分校
有人说 R 不适合高性能计算 (HPC)。 也有人说：“不……等等……你是在开玩笑吗？R 用于 HPC？”
然而，世界正在改变。数据分析是新的热门、酷炫的东西。无论您是看到数据科学带来的切实利益，还是仅仅看到利润，事实是 R 在这方面表现出色。HPC 的格局正在发生变化，以更好地适应数据分析应用程序和用户。因此，由于 R 压倒性的受欢迎程度，它是您应该关注的自然选择。幸运的是，许多人已经响应号召，将 R 嵌入到 HPC 环境中。
我在这里对我的听众做了一些假设。我猜您可能对 R 了解不多，但至少感到好奇。也许这种好奇心是由应用程序需求驱动的，害怕错过（所有酷炫的人都在用 R 编程），或者您可能有越来越多的客户向您咨询 R 解决方案。无论您的动机是什么，欢迎您。我很高兴有您。我还假设您在 HPC 方面相当精通。您知道什么是编译器，并且您最喜欢 Intel® 编译器。
基于这些假设，本文将有所不同。我想向您介绍基础知识，并介绍 R 领域中的一些竞争性想法，同时尽量不偏不倚。也许谈论在数千个节点上使用 R 进行分布式计算，并获得 TB 数据的交互式速度更有趣——我们可以做到。但 R 在 HPC 界相对不为人知，我感到有必要在“像炮弹一样冲出去”之前花点时间学习“爬行”。
我们将从一点历史和一些基础知识开始，这些知识是任何接触 R 并关心性能的人都应该了解的。我们将花更多时间讨论将编译代码集成到 R 中，然后讨论并行计算。这有点像一阵旋风，本文不会让您成为任何一个主题的专家。但希望能够为您认真开始使用 R 提供足够的帮助。
背景
R 的根源可以追溯到 1976 年，当时贝尔实验室的 John Chambers 开始研究 S。S 最初被设计为一系列 Fortran 代码的交互式接口。而且，无论您多么努力，您都无法摆脱 Fortran，所以 R 今天也基本是这样工作的。R 本身于 20 世纪 90 年代初发布，由 Ross Ihaka 和 Robert Gentleman 创建，是一个免费的 S。严格来说，R 是 S 语言的一个方言，这是一个非常自命不凡的说法，意味着 80 年代编写的很多 S 代码仍然可以在 R 中运行。
R 是一种奇怪的语言。我最喜欢用来展示这一点的一个例子是
typeof(1) ## [1] "double" typeof(2) ## [1] "double" # 1:2 is the vector of numbers "1 2" typeof(1:2) ## [1] "integer"
如果您盯着它看足够长的时间，它似乎有点道理。很自然地可以假设 1:2 可能是一个索引。但它仍然很奇怪。通常，: 有点不可预测。如果您执行“1”:“2”，它将返回整数向量 1 和 2。所以它只是对字符进行 ASCII 转换并转换为整数，对吧？但“A”:“B”会报错。而且 : 并不总是产生整数；例如，您可以执行 1.5:2.5。
R 拥有一个非常活跃的社区，其综合 R 存档网络 (CRAN) 上拥有超过 10,000 个贡献的包。这些包涵盖了从硬核数值计算，到尖端的统计和数据科学方法，再到构建交互式数据分析网站（称为Shiny，它非常棒）。就我而言，R 的包基础设施是无与伦比的。它几乎是您能找到的“它就是能用”的东西。所以，对于所有那些一直窃笑并点头称“确实如此”的 Python* 粉丝来说，现在可能是时候解释为什么一群您认为不知所云的统计学家竟然设法创建了唯一一个不是糟糕透顶的打包框架了。
现在，当我 S R 很受欢迎时，我没有开玩笑。2016 年，IEEE Spectrum 的编程语言排名将 R 排在第五位，超过了 C# 和 JavaScript*。尤其有趣的是，这是对编程语言的排名。即使是那些热爱 R 的人也会告诉你，它是一种糟糕的编程语言。R 在数据分析方面非常出色，以至于人们愿意忽略它所有的怪癖，去发现隐藏在下面的真正美丽之处。
换句话说，R 就像《加勒比海盗》电影中的杰克·斯派洛：它可能是您听说过的最差（语言）……但您听说过它。
免费改进
美国喜剧演员 W. C. Fields 曾说过：“我遇到的最懒惰的人，把爆米花放进煎饼里，这样它们就会自己翻面。” 我想 Fields 所说的这个传说中的人会成为一名出色的工程师。毕竟，既然别人可以替我们做繁重的工作，又何必辛苦呢？
在 R 中，有几种方法可以进行这种“站在巨人的肩膀上”。首先——这应该并不令人意外——如果您使用好的编译器编译 R，您可以期望看到一些显著的性能提升。R 是用 C、Fortran 和 R 编写的，因此在 Intel® 硬件上使用 Intel 的 icc 和 ifort 是一个不错的起点。而且，幸运的是，Intel 有一篇非常好的文章，介绍了如何在 Linux 上使用 Intel 编译器构建 R。
这是从 R 中获得良好性能的有力第一步。但 R 中构成基础 R 的所有 R 代码怎么样？R 自 2.13.0 版本以来就有一个字节码编译器，并且它为 2.14.0 及更高版本编译了 R 内部代码。对于您自己编写的代码，传统上您需要付出一些努力才能使用字节码编译器。然而，在 3.4.0 版本（在本篇撰写之时即将发布）中，R 将包含一个 JIT，这将使许多关于使用编译器的旧建议变得不那么重要。
现在，值得指出的是，字节码编译器并不像真正的编译器那么好。如果您的代码设计不当，例如不必要地进行计算，那么它仍然在那里；计算只是以字节码形式执行。它也不是 R 到 C 的翻译器之类的东西。它在循环密集型代码（不包括隐式循环）上效果最好，（从性能角度来看）否则几乎什么作用都没有。我曾见过它使循环体性能提高约 10%，也见过它对性能的影响只有 0.01%。但嘿，这不需要您付出任何努力，所以我们姑且接受。
这些改进都很好，肯定有助于提高 R 代码的运行时间，但不会让您欣喜若狂。现在，如果您想买一对新袜子，那么选择好的 LAPACK 和 BLAS 库可以带来非常令人印象深刻的性能提升。这些是矩阵运算的事实标准数值库，R 使用它们来驱动其低级线性代数以及大部分统计运算。讽刺的是，统计学中最重要的运算之一，线性回归，并不使用 LAPACK。相反，它使用 LINPACK 的一个高度修改的版本。不，不是运行在超级计算机上的基准测试。我说的是 LAPACK 的 70 年代的前身。原因有点复杂，但确实存在原因。所以您经过优化的 LAPACK 对线性回归没有帮助，但仍然可以利用好的 Level-1 BLAS。
R 附带所谓的“参考”BLAS，它速度慢得令人发指。尽管有上述例子，但如果您将 R 与好的 BLAS 和 LAPACK 实现链接，那么您可以期望看到显著的性能改进。而且，碰巧的是，Intel 在Intel® 数学核心库 (Intel® MKL) 中有一个非常高质量的实现。据我所知，Microsoft 提供了一个 R 版本，它是用 Intel 编译器编译的 R，并附带 Intel MKL。他们称之为Microsoft R Open，它是免费提供的。他们还维护了一个详细的基准测试集合，展示了 Intel MKL 的威力。或者，如果您愿意，您可以遵循 Intel 自己的文档，了解如何将 R 与 Intel MKL 链接。
对于那些拥有所有最新、最炫的玩具的人：是的，这也适用于 MIC 加速器。其中许多早期工作来自德克萨斯高级计算中心 (TACC) 的优秀人士，他们对使用 Intel MKL 自动卸载进行了大量实验。说事情进展顺利有点轻描淡写了，拥有少量 Intel® Xeon Phi™ 处理器的 R 用户应该认真考虑尝试一下。如果您不确定从哪里开始，Intel 还提供了一份非常实用的指南，帮助您处理这类事情。
利用编译代码
R 领域正在发生的有趣革命之一是 C++ 在 R 包中的使用日益增加。这主要归功于 Dirk Eddelbuettel 和 Romain Francois，他们创建了Rcpp 包。Rcpp 极大地简化了将 C++ 代码集成到 R 分析管道中的过程。是的，他们设法说服了一群认为 Python 太复杂而无法用 C++ 编程的统计学家。我和您一样惊讶。但无论他们如何做到这一点，我们都是受益者。这意味着 CRAN 包的速度越来越快，内存消耗越来越少。而使用 Intel 编译器的人将从这场革命中受益最多，因为正如我们之前讨论过的：更好的编译器，更快的代码。
现在，对于那些喜欢标准 C 的人，R 提供了一个一流的 C API。事实上，Rcpp 就是在此基础上构建的；尽管我认为 Rcpp 付出了更大的努力是公平的。但回到 C API，对于那些能够容忍 Fortran 并愿意编写 C 包装器而不希望 C++ 链接器参与其中的人来说，这也很有用。该 API 主要记录在编写 R 扩展手册中，这是任何从事 R 工作的人不可或缺的资源。但是，为了回答一些出现的问题，如果您走这条路，您可能会发现自己需要查看 R 的头文件。
对于 Rcpp 路线，您像安装其他任何 R 包一样安装它，即
install.packages("Rcpp")
为了快速了解这个包的威力，让我们看看“数值 Hello World”，即您已经见过无数次的用于查找 Π 的蒙特卡罗积分。在 R 中，您可能会这样写：mcpi <- function(n)
mcpi <- function(n)
{
r <- 0L
for (i in 1:n){
u <- runif(1)
v <- runif(1)
if (u*u + v*v <= 1)
r <- r + 1L
}
return(4*r/n)
}
现在，当我 S“您可能会写”时，我再次假设您对 R 并不熟悉。可能没有任何有经验的 R 用户会写这样的东西。人们可以合理地争辩说，它看起来有点像 C。我给任何将要维护 R 代码库以提高其速度的人最好的建议是：它越像 C，在 R 中运行得越差，但越容易转换为 C/C++。反之亦然，它越不像 R，转换起来就越困难。更自然的 R 解决方案是以下向量化乱码
mcpi_vec <- function(n)
{
x <- matrix(runif(n * 2), ncol=2)
r <- sum(rowSums(x^2) <= 1)
return(4*r/n)
}
与所有其他高级语言一样，使用向量化https://software.intel.com/en-us/intel-vectorization-tools>将提高运行时性能，但也会消耗更多内存。所以，让我们忘记所有这些 R 的事情，只写一个 C++ 版本
#include <Rcpp.h>
// [[Rcpp::export]]
double mcpi_rcpp(const int n)
{
int r = 0;
for (int i=0; i<n; i++){
double u = R::runif(0, 1);
double v = R::runif(0, 1);
if (u*u + v*v <= 1)
r++;
}
return (double) 4.*r/n;
}
现在，除了那个神秘的 Rcpp::export 部分，它看起来可能非常易读的 C++ 代码。事实证明，神秘的部分将处理所有“样板”代码的生成。在 R 中，没有标量（嘿，我告诉过您它是一种奇怪的语言），只有长度为 1 的向量。因此，在后台，Rcpp 实际上会为您处理这种心智负担，并在其包装器中为您创建一个长度为 1 的双精度向量。通常，我们可以使用整数和双精度数以及这些基本类型的向量来玩这个游戏。更复杂的事情涉及更多的复杂性。但嘿，这不是很棒吗？
要编译/链接/加载并生成各种样板代码，我们只需要调用 sourceCpp() 即可使该函数立即在 R 中可用
Rcpp::sourceCpp(file="mcpi.cpp")
mcpi_rcpp(10000)
## [1] 3.1456
眼尖的读者可能会问：“这里的‘10000’难道不是双精度数吗？”您想得没错，因为它确实是。我们可以通过调用 10000L 来要求一个整数——请注意，这是一个普通的 32 位整数——但 Rcpp 会自动为您处理类型转换。它实际上处理的转换与 R 代码版本完全相同。这有利有弊，但这是他们采用的方法，值得注意并了解。
我们可以使用 rbenchmark 或 microbenchmark 包轻松比较这三者的性能。我个人比较喜欢 rbenchmark，它大致是这样的
library(rbenchmark)
n <- 100000
cols <- c("test", "replications", "elapsed", "relative")
benchmark(mcpi(n), mcpi_vec(n), mcpi_rcpp(n), columns=cols)
## test replications elapsed relative
## 1 mcpi(n) 100 49.901 214.167
## 2 mcpi_vec(n) 100 1.307 5.609
## 3 mcpi_rcpp(n) 100 0.233 1.000
而且，这很不错！现在，当然，这为 OpenMP* 或Intel® 线程构建块 (Intel® TBB) 等提供了机会。但说到并行……
并行编程
自 2.14.0 版本以来，R 就附带了 parallel 包。它通过提供两个独立的 API（一个使用套接字，一个使用 OS fork）来实现非常简单的任务级并行。存在两个接口一部分是历史原因，因为它们源自较旧的贡献包 multicore 和 snow。但保留两者的愿望最好通过 R 核心希望支持所有平台，即使是缺乏 fork 的 Windows* 来解释。在非 Windows 平台上，感兴趣的函数是 mclapply()，它是 multicore lapply()——之所以这样命名，是因为它应用一个函数并返回一个列表。在这里，R 展示了它的一些函数式编程能力
lapply(my_data, my_function)
parallel::mclapply(my_data, my_function)
数据可以是索引，也可以是包含非常大的、复杂对象的卷积列表。只要提供的函数可以处理输入，它就会起作用。
这是两个官方支持的接口之一。另一个更复杂，通常只供 Windows 程序员使用。这在 R 用户之间造成了一些分歧。这纯粹是我的个人观点，但我认为 R 用户不像其他语言中的普通程序员那样喜欢有多种选择。他们想要一种好的做事方式，然后就结束了。为此（讽刺的是），已经出现了几个项目来尝试统一所有分散的接口。其中包括较旧且更成熟的 foreach 包，以及来自 Bioconductor 项目的较新的 BiocParallel。
您可能会想，为什么有两个独立的接口很重要。实际上，有更多的包可以在 R 中实现并行。HPC 任务视图https://cran.r-project.cn/web/views/HighPerformanceComputing.html 是发现众多选项的好资源。
如果您觉得大部分关注和兴趣都集中在共享内存并行上，那么您就对了。总体而言，R 社区在相对较近的时期之前一直有点抗拒关心性能。我认为这很大程度上是因为 R 的思维空间仍然被数据科学的统计学方面所主导。坦率地说，统计学中的大多数大数据问题并非如此。您仍然可以通过对数据进行降采样并使用经典的统计技术来完成很多工作。它不是万能的，但肯定有其用武之地——而且这些工具可能被低估了。
但这都不涉及超级计算机，所以忘掉那些胡说八道吧。让我们谈谈 MPI。Rmpi 包可以追溯到 2002 年。最近，R 中的大数据编程 (pbdR) 项目一直在开发用于在超级计算环境中使用 R 进行大规模计算的包。现在，坦白说：我参与了这个项目。因此我的观点自然带有偏见。但我认为我们有一些有趣的东西可以展示给您。
我们维护了相当多的包，并且通常致力于将 HPC 的最佳实践带入 R，以用于数据分析和性能分析。为简洁起见，让我们只关注直接的 MPI 编程。我们将简要比较 Rmpi 和我们的包 pbdMPI。首先，也是最重要的一点，Rmpi 可以交互使用，但 pbdMPI 不能。这是因为 pbdMPI 被设计为仅用于单程序多数据 (SPMD) 风格，而 Rmpi 实际上是用于经理/工作者风格。如果您曾经在集群上使用 MPI 提交批处理作业，那么您几乎肯定是在编写 SPMD。这是一个如此直观的想法，如果您以前从未听说过，您会惊讶于它竟然有一个名字。因此，对于那些从 HPC 世界来到 R 的人来说，pbdMPI 应该感觉很熟悉。
除了编程风格，两个包的 API 之间也存在一些严重差异。在 Rmpi 中，您需要指定数据的类型。例如
library(Rmpi)
mpi.allreduce(x, type=1) # int
mpi.allreduce(x, type=2) # double
还记得我们开头处的类型示例吗？在 pbdMPI 中，我们使用 R 的面向对象功能来尝试自动处理这些以及其他底层细节
library(pbdMPI)
allreduce(x)
这是一个小例子，但很好地展示了我们的理念。我们认为 HPC 社区做得非常出色，但我们也认为 HPC 工具对大多数人来说太难使用了，应该变得更简单。
为了一个稍微实质性的例子，让我们快速看一下并行“Hello World”。我们提到 pbdMPI 必须在批处理模式下使用。缺点是 R 用户很难将批处理而不是交互式处理作为一种思维方式。优点是它能很好地与您已经知道的所有 HPC 功能（如资源管理器和作业调度程序）协同工作。假设我们想运行我们的“Hello World”示例
library(pbdMPI)
comm.print(paste("Hello from rank", comm.rank(), "of", comm.size()), all.rank=TRUE)
finalize()
我们只需要调用适当的 mpirun（或您系统的等效命令）
mpirun -np 2 Rscript hello_world.r
这将产生预期的输出
[1] "Hello from rank 0 of 2"
[1] "Hello from rank 1 of 2"
摘要
统计学界有一个关于 George Box 的名言：所有模型都是错误的，但有些是有用的。那么，我认为所有编程语言都是糟糕的，但有些是有用的。R 可能是这个想法的终极体现。毕竟，如果它被坚持了 40 年（算上 S），并且目前在编程语言中排名第五，那一定有它的道理。虽然 R 以速度慢而闻名，但肯定有一些策略可以缓解这种情况。使用好的编译器。好的 BLAS 和 LAPACK 将提高许多数据科学操作的性能。将编译后的内核嵌入到您的 R 分析管道中可以极大地提高性能。而且，在不确定的时候，就用更多的核心来解决您的问题。