1. 问题背景

在目前的数据科学（data science）领域中，R是一个被广泛采用的统计计算语言。R因其丰富的统计计算功能、友好的交互式用户界面和易于使用的线性代数、数据分析API而深受数据分析师的喜爱。但R语言及其计算环境最初是为单机串行计算设计的，R语言解释器采用了全单机内存的处理方式，这种设计限制R可处理的数据集大小。另外，R语言默认的是单线程实现，无法利用多核计算算能力。为了使R语言能够处理大数据，一些工作开始尝试将R语言计算环境与大数据计算系统Hadoop、Spark相融合。但遗憾的是这些方案都缺乏对线性代数计算的支持，并且他们采用了不同于R本身的编程模型，用户需要改写已有的R程序才能完成并行计算。目前在R语言的软件生态中，缺乏相关的软件包或系统，使用户可以交互式地进行大规模分布式线性代数与数据分析计算。

本项目提出了一个系统iPLAR，用以弥补交互式R语言环境中无法进行大规模分布式线性代数与数据分析计算的不足。iPLAR借助Client/Server架构，将交互式R语言环境与基于MPI的并行计算环境相解耦，使用户可以在交互式的R语言环境中，利用MPI提供的分布式计算能力来加速计算。在MPI端，iPLAR采用ScaLAPACK库作为高性能分布式线性代数计算软件库。在面向用户的交互式R环境中，iPLAR对底层实现进行了封装，向用户提供了一套与R原生线性代数API类似的用户接口，用户能以较小的移植代价将原生的R程序移植为iPLAR程序，而不需要了解MPI和ScaLAPACK的编程细节，达到了MPI编程对用户的透明性。

2. 问题分析

本工作的主要的核心问题是解决交互式R语言环境与批处理式MPI线性代数计算环境之间的解耦，并在两个计算环境之间搭建良好的交互机制。

3. 系统总体架构/算法总体流程

3.1 iPLAR的软件架构

iPLAR系统是基于一系列R软件生态中的其他软件建立的。iPLAR的软件架构如图1所示。Client端运行在交互式R计算环境中，向用户提供一组用于操作分布式稠密矩阵的API。Client端采用面向对象技术，将对用户提供的API都包装在iddmatrix（interactive distributed dense matrix）类型的对象中。用户面向iddmatrix对象进行编程。每个iddmatrix对象都代表着一个远程的分布式稠密矩阵。iddmatrix重载了原生R矩阵的方法，使其API与R原生矩阵兼容。而在Server端，Server端最上层是一个与Client端进行交互的Server Interface。Server Interface负责解析计算指令并执行指令。Server Interface在执行过程中具体的会调用pbdDMAT软件包进行分布式矩阵计算。pbdDMAT是pbdR软件项目中提供的一个软件包，为R语言用户提供了对ScaLAPACK软件库的一层抽象封装。pbdDMAT本身是基于MPI的SIMD编程模型开发的，只能在MPI计算环境中使用，无法直接应用于交互式R语言会话。iPLAR通过对pbdDMAT的封装和间接调用，使其提供的计算方法能在交互式R环境中使用。pbdDMAT在底层是通过调用ScaLAPACK软件库来实现计算，而ScaLAPACK在下层依次调用PBLAS、BLACS等软件库来完成细粒度的线性代数计算任务。Server端通过HDFS来提供分布式矩阵的文件存储。

3.2 iPLAR的系统架构

iPLAR采用Client/Server架构，以将交互式的R语言环境与批处理式的MPI计算环境进行隔离。iPLAR的模块架构和运行流程如图1所示。Client端运行在交互式R语言环境中，负责与R用户进行交互；Server端负责分布式地矩阵存储与计算。客户端和服务器之间通过R Raw Socket机制完成网络通信。除了Client、Server两个模块，iPLAR系统还包含有一个Daemon模块。Daemon模块负责与集群作业管理器交互来完成MPI计算进程的启动。

从用户的角度看，R用户只需在其本地的交互式R环境中加载Client端的软件包，便可以像使用原生R矩阵一样对分布式矩阵进行线性代数、数据分析类运算，Client端内部与Server端的交互过程对用户完全透明，用户无需了解任何分布式编程的知识即可使用。

图2 iPLAR系统架构与运行流程图

4. 技术实现方案与优化

4.1 Client端设计

Client端主要负责在用户与iPLAR系统之间架起交互的桥梁。iPLAR的Client端设计主要采用面向对象编程技术。其提供了一个R语言的S4类iddmatrix类，作为面向用户的API。每个iddmatrix类的对象都代表了Server端的一个分布式稠密矩阵对象（对应于ddmatrix对象）的操作句柄。在Client端只保存远程分布式矩阵的名称和尺寸信息，而不保存具体的矩阵数据。所有对Client端iddmatrix对象的操作，都会转换为计算指令发送给Server端；Server端在解析指令并执行后，将结果矩阵的句柄返回给Client端；Client端在拿到句柄后封装为一个新的iddmatrix对象返回给用户。Client端与Server端之间只传递操作句柄，并不传递具体数据；具体的矩阵数据在Server端以ScaLAPACK矩阵的形式分布式地存储在内存中。通过这样的设计，Client端避免了存储大量的具体矩阵数据，使Client端能操作远超本地内存大小的超大规模矩阵。

Client端重载了iddmatrix类的线性代数计算操作符，使iddmatrix类对象的方法接口与R语言原生提供的矩阵对象的方法接口保持一致。Client端实现的iddmatrix类的操作符列表见表1，函数的命名与参数规范与R语言原生的函数定义保持一致。这样iPLAR系统与R原生矩阵在绝大部分的操作上保持了API级兼容。Client端面向对象的封装使Client端与Server端交互的细节对用户完全透明，除了启动iPLAR会话的操作之外，用户感知不到Server端的存在，实现了用户与ScaLAPACK/MPI编程的隔离，使ScaLAPACK库能以一种间接的方式在交互式的R语言环境中被使用。而iddmatrix与原生矩阵API级的兼容使用户的学习成本很小，将已有的R程序迁移到iPLAR系统上的迁移代价很小。图4展示了使用iPLAR编写Logistic Regression训练算法的编程示例。可以发现在27行的代码中，iPLAR版的代码与原生R代码相比只有1行不同，移植开销很小。另外在BP神经网络训练算法的实现中 ，60行原生R代码只需要修改3行即可移植为iPLAR版的代码。上述示例说明iPLAR系统学习代价小，移植开销小。

4.2 Server端设计

iPLAR系统的Server端是一个面向大规模矩阵计算的分布式计算引擎。它负责矩阵数据的具体存储、线性代数计算和数据并行计算。它支持分布式矩阵的按行、按列、按块三种存储方式，以满足不用计算需求。Server端的执行流程类似“读取-求值-输出”（REPL）的循环：

1. 等待Client端发送计算指令；
2. 接收到指令后，解析Client端计算指令，并执行计算；
3. 如果计算过程生成新的分布式矩阵对象，为其分配句柄；
4. 将计算结果信息（成功、失败）和结果分布式矩阵的句柄返回给客户端；
5. 跳转第1步执行。

Client端与Server端交互所使用的指令格式如图3所示。一条指令对应于一个R语言的列表（list），通过R Raw Socket进行传递。列表的GList域指定了垃圾矩阵对象的句柄列表，用于Client端与Server端之间的同步垃圾回收；REQ.TYPE域指定了计算指令的类型，Server端据此决定该执行何种操作；Body域保存了指令的操作数信息，即参与运算的分布式矩阵句柄，Server端根据该句柄，确定需要参与运算的分布式矩阵对象，然后调用pbdDMAT提供的对应API完成基于ScaLAPACK库的分布式线性代数计算，或者使用数据并行API完成数据分析计算。

5. 性能评估分析

本节所述实验是在一个有10个计算节点组成的集群中进行的。每个计算节点的计算性能如表2所示。实验过程中将使用ATLAS[9]作为基础的BLAS运算库的实现。所有针对单机R环境的测试中，均使用多线程版本的ATLAS库以充分利用单机的性能；所有针对iPLAR的计算测试中均采用多进程并发的方式，每个进程中采用单线程版本的ATLAS库，每台机器最多启动8个进程。

在本节所展示的实验中，iPLAR的性能将与R原生提供的单机计算能力、基于核外计算技术的big.matrix软件包提供的核外线性代数计算能力进行对比。Big.matrix的实现来自于CRAN上的“Big Memory”系列软件包 ，该软件包提供了基于核外计算的线性代数API，使用户可以借助磁盘存储来处理比单机内存更大的矩阵数据。Big.matrix、原生R以及iPLAR都提供了交互式环境下的线性代数计算能力，其他的方案要么需要大规模重写用户程序（例如使用Ricardo或SparkR）、要么不支持交互式的使用（例如pbdDMAT）。

5.1 基本线性代数操作性能测试

根据不同的线性代数函数计算过程中纯计算与纯通信开销的比例不同，选取3种典型的线性代数基本操作用于性能测试与评估：

1. 计算密集型操作：按元素求自然指数（exp）；
2. 通信密集型操作：矩阵转置（t）；
3. 计算通信双密集型操作：矩阵乘法（%*%）。

为了测试iPLAR的C/S架构引入的额外开销，我们测试了一个空操作所耗的时间。在一次空操作中，除了没有任何计算过程之外，其他流程与普通的计算操作完全相同。在MPI的Rank数从16变到80的过程中，空操作所耗时间从0.07s涨到0.11s。一个空操作所耗时间可以认为是iPLAR引入的整套架构所带来的额外开销，与后面测试中可以看到的实际计算开销相比，几乎可忽略不计。使用iPLAR进行计算时，主要时间开销还是花费在具体的ScaLAPACK计算上。

对于iPLAR的数据可扩展性性能，实验比较了iPLAR、单机多线程原生R、基于核外计算技术的big.matrix软件包提供的线性代数计算能力，计算结果如图4所示。可以发现iPLAR具有近似线性的数据可扩展性，而且性能比基于核外计算技术的big.matrix要好。对于非纯通信密集型的计算（exp和矩阵乘法）中，iPLAR获得了比单机多核原生R、以及big.matrix更好的计算性能。对于通信密集型计算（矩阵转置），iPLAR的分布式实现受限于网络通信速度，性能落后于单机的实现，但仍然比基于核外计算的big.matrix要好。

对机器可扩展性的评估再一次印证了数据可扩展性实验中表现出的结果：对于非通信密集型的计算，iPLAR提供了良好的机器可扩展性；但对于纯通信密集型的计算，则受限于网络通信，节点数越多计算开销越大。具体结果如图5所示。

5.2 机器学习算法性能测试

本节选取两个典型的机器学习分析算法作为综合测试用例，来综合考察iPLAR在实际应用中所体现的性能。选取的两个应用是基于梯度下降的Logistic Regression训练算法（LR）和BP神经网络训练算法（BPNN），这两个训练算法分别是线性模型和非线性模型训练的典型代表。

图6展示了两个综合应用的数据可扩展性实验结果。因为big.matrix无法完成按元素自然对数运算（exp），故其无法完成BPNN算法，所以实验结果缺失。iPLAR展现出了近似线性的数据可扩展性；借助iPLAR，在10节点的集群中，用户可以交互式地处理比单机计算极限大10倍的数据集。对于轻量型的训练算法（LR），因为其主要操作多为计算密集型操作，所以他取得了比单机R和big.matrix更好的性能。但是对于通信密集型操作较多的算法（BPNN）iPLAR的性能就不如单机R的性能，但是借助分布式计算，iPLAR可以处理比单机内存极限大的多的数据集。

图7展示了机器可扩展性的实验结果，在两个综合测试用例中，iPLAR的机器可扩展性均表现一般。对于6个节点时的性能变差的原因主要是处理器网格的自动设置不当，导致性能变差。ScaLAPACK里的处理器网格设置对计算性能有明显影响，而如何自动化地选择良好的设置仍是一个开放的研究问题。综合上述表现，iPLAR对于数据科学家和分析师的意义将会主要在于数据扩展性上，用户可以利用iPLAR提升交互式环境的计算能力、无缝的处理更多的数据而无需感知分布式计算的存在；在某些算法（例如LR）中，分布式并行处理能带来更快的性能。

图 7 综合应用的机器可扩展性实验结果。训练样本集大小均为\(10^7\)。

6. 总结

iPLAR在交互式R语言计算环境与基于MPI的分布式线性代数计算软件ScaLAPACK之间架起了一座桥梁，使普通的R语言用户可以在不需要了解MPI和ScaLAPACK编程细节的情况下，就能利用起分布式计算的能力，快速高效地完成线性代数计算任务和基于线性代数的数据分析任务。iPLAR采用了Client/Server模式的软件架构。在Client端，iPLAR为用户提供了基于线性代数计算的向量化编程模型，并且API与R原生的线性代数API保持一致，这样保持了较小的学习代价和迁移代价。对iPLAR的性能评估表明，借助分布式并行计算技术，在一个10节点的集群中，iPLAR能处理的最大数据量是单机R环境的10倍；iPLAR提供了近线性的数据可扩展能力。为了克服全局向量化编程在大数据分析中的潜在缺陷，iPLAR同时也提供了将局部向量化编程与全局数据并行相结合的混合数据并行模型。通过混合数据并行编程，两个典型的机器学习算法（Logistic 回归与BP神经网络训练）能获得近线性的数据可扩展性和良好的机器可扩展性。

7. 工作产出

1. 论文：Zhaokang Wang, Shiqing Fan, Rong Gu, Chunfeng Yuan, and Yihua Huang, "iPLAR: Towards Interactive Programming with Parallel Linear Algebra in R." In Algorithms and Architectures for Parallel Processing, vol. 9531, pp. 104–117. Springer LNCS, 2015. (ICA3PP 2015)
2. 专利：在交互式R语言平台中进行并行线性代数计算的方法。申请号: 201510755923.2，申请公布号: 申请公布号：CN105389220A。