Large Scale Learning in R

Wush Wu
Taiwan R User Group

Taiwan R User Group

MLDMMonday: 每周一分享资料相关议题

主题包含但不限于:

R 套件使用

机器学习和统计模型

很多人都在诟病R 无法处理大量的数据

但是只要用对R 的包

R 是可以处理大量的数据

今天我想跟大家分享

运用R 建立商用推荐引擎的心得

主要是靠着RcpppbdMPI等包

打造可扩放性的学习模组

实际的数据和课本上的数据是不一样的

实际的数据

  • 乱、充满错误
  • 不停的变动
  • 不知道怎样才算好
    • 只能不断精益求精

课本的数据

  • 干净
  • 静止
  • 需要复杂的算法

今天从这样的数据开始

is_click show_time client_ip adid
1 FALSE 2014/05/17 04:06:52 114.44.x.x 133594
2 FALSE 2014/05/17 04:04:48 1.162.x.x 134811
3 FALSE 2014/05/17 04:03:05 101.13.x.x 131990
4 FALSE 2014/05/17 04:05:04 24.84.x.x 134689
5 FALSE 2014/05/17 04:03:26 140.109.x.x 126437
6 FALSE 2014/05/17 04:04:28 61.231.x.x 131389


以“预测点击发生的机率”为例

问题的建模

is_click show_time client_ip adid
1 FALSE 2014/05/17 04:06:52 114.44.x.x 133594
2 FALSE 2014/05/17 04:04:48 1.162.x.x 134811
3 FALSE 2014/05/17 04:03:05 101.13.x.x 131990
4 FALSE 2014/05/17 04:05:04 24.84.x.x 134689
5 FALSE 2014/05/17 04:03:26 140.109.x.x 126437
6 FALSE 2014/05/17 04:04:28 61.231.x.x 131389

$$\Downarrow$$

$$P(y|x) = \frac{1}{1 + e^{-yw^Tx}}$$

\[w^Tx\]

iris

head(model.matrix(Species ~ ., iris))
  (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
1           1          5.1         3.5          1.4         0.2
2           1          4.9         3.0          1.4         0.2
3           1          4.7         3.2          1.3         0.2
4           1          4.6         3.1          1.5         0.2
5           1          5.0         3.6          1.4         0.2
6           1          5.4         3.9          1.7         0.4

实际的资料

有 \(10^{9+}\) 笔资料

有 \(10^{3+}\) 笔广告

每种类别变量可能到 \(10^{3}\) 类别

model.matrix的后果:

至少需要 \(10^9 \times (10^3...)\) 个位元

也就是 \(116.415\) GB的记忆体

开始抱怨了!

R 不能处理大量的数据

R 会吃很多记忆体

使用R 跑不动,使用其他工具就跑得动

请使用和其他工具一样的资料结构:

这种状况,应该使用稀疏矩阵

稀疏矩阵

  • 只储存非0的资讯
n <- 1
m1 <- matrix(0, 10^n, 10^n);m1[1, 4] <- 1
c(m1)
library(Matrix)
m2 <- Matrix(0, 10^n, 10^n, sparse=TRUE)
m2[1,4] <- 1
str(m2)
  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [48] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [95] 0 0 0 0 0 0
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int 0
  ..@ p       : int [1:11] 0 0 0 0 1 1 1 1 1 1 ...
  ..@ Dim     : int [1:2] 10 10
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num 1
  ..@ factors : list()

在我手上的资料上

大概可以省下 \(10^2 \sim 10^3\) 倍记忆体

而且可以大幅度加快运算效能!

如果m1, m2是 \(10^4 \times 10^4\) 的矩阵:

      test replications elapsed relative user.self sys.self user.child sys.child
1 m1 %*% r          100  21.690      868    47.727    0.217          0         0
2 m2 %*% r          100   0.025        1     0.021    0.003          0         0

可以利用Rcpp将C的实作暴露到R中

可以利用Rcpp做记忆体的重复使用

                test replications elapsed relative user.self sys.self user.child sys.child
1           m2 %*% r          100   0.023     7.67     0.021    0.002          0         0
2 XTv(m2, r, retval)          100   0.003     1.00     0.002    0.000          0         0

在现代的Rcpp架构下

将C++函数放到R中变得很简单

我们只需要专注在算法上

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
SEXP XTv(S4 m, NumericVector v, NumericVector& retval) {
  //...
}

再透过pbdMPI开发分散式矩阵乘法

利用更多CPU和更多记忆体提升效能

\[\left(\begin{array}{c}X_1 \\ X_2\end{array}\right) v = \left(\begin{array}{c}X_1 v \\ X_2 v\end{array}\right)\]

\[\left(v_1 , v_2\right) \left(\begin{array}{c}X_1 \\ X_2\end{array}\right) = v_1 X_1 + v_2 X_2\]

MPI

记忆体足够的话,较快

没有容错

Hadoop

慢,要大量机器才有效果

有容错

pbdMPI的优点

好安装

好开发

library(pbdMPI)
.rank <- comm.rank()
filename <- sprintf("%d.csv", .rank)
data <- read.csv(filename)
target <- reduce(sum(data$value), op="sum")
finalize()

最佳化算法

迭代次数少

不用计算Hessian矩阵

已有很棒的实作LIBLINEAR

为了加强效能

并且能够更方便的更改模型

我们自己包LIBLINEAR到R中

LIBLINEAR 原始码中的 tron.h

class TRON
{
public:
  TRON(const function *fun_obj, double eps = 0.1, int max_iter = 1000);
  ~TRON();

  void tron(double *w);
    void set_print_string(void (*i_print) (const char *buf));

private:
    //...
};
  • TRON是最佳化的核心算法的实作
  • TRON没有牵涉到资料,甚至没有牵涉到Loss Function
  • function提供了一个界面

界面:function

class function
{
public:
  virtual double fun(double *w) = 0 ;
  virtual void grad(double *w, double *g) = 0 ;
    virtual void Hv(double *s, double *Hs) = 0 ;

    virtual int get_nr_variable(void) = 0 ;
    virtual ~function(void){}
};
  • fun代表objective function
  • gradfun的gradient
  • Hvfun的hessian乘上一个向量

利用Rcpp Modulesfunction暴露到R中

实作Rfunction

class Rfunction : public ::function {
  Rcpp::Function _fun, _grad, _Hv;
  Rcpp::NumericVector _w, _s;

  Rfunction(SEXP fun, SEXP grad, SEXP Hv, int n) 
  : _fun(fun), _grad(grad), _Hv(Hv), _w(n), _s(n) { }

  //...
};

SEXP tron//...

将Rfunction暴露出来

RCPP_MODULE(HsTrust) {  
  class_<Rfunction>("HsTrust")
    .constructor<SEXP, SEXP, SEXP, int>()
    .property("n", &Rfunction::get_nr_variable, 
    "Number of parameters")
    .method("tron", &tron)
    .method("tron_with_begin", 
    &tron_with_begin)
    ;   
}

实作的结果可包装成套件

library(devtools)
install_bitbucket(
  repo="largescalelogisticregression", 
  username="wush978", ref="hstrust")

使用

library(HsTrust)

beta <- 1
fun <- function(w) 
  sum((w-beta)^4)

grad <- function(w)
  4 * (w-beta)^3

Hs <- function(w, s) 
  12 * (w-beta)^2 * s

obj <- new(HsTrust, fun, grad, Hs, 1)
print(r <- obj$tron(1e-4, TRUE))
library(HsTrust)
# ...
print(r <- obj$tron(1e-4, TRUE))
iter  1 act 8.025e-01 pre 6.667e-01 delta 4.186e-01 f 1.000e+00 |g| 4.000e+00 CG   1
iter  2 act 1.585e-01 pre 1.317e-01 delta 4.186e-01 f 1.975e-01 |g| 1.185e+00 CG   1
iter  3 act 3.131e-02 pre 2.601e-02 delta 4.186e-01 f 3.902e-02 |g| 3.512e-01 CG   1
iter  4 act 6.185e-03 pre 5.138e-03 delta 4.186e-01 f 7.707e-03 |g| 1.040e-01 CG   1
iter  5 act 1.222e-03 pre 1.015e-03 delta 4.186e-01 f 1.522e-03 |g| 3.083e-02 CG   1
iter  6 act 2.413e-04 pre 2.005e-04 delta 4.186e-01 f 3.007e-04 |g| 9.135e-03 CG   1
iter  7 act 4.767e-05 pre 3.960e-05 delta 4.186e-01 f 5.940e-05 |g| 2.707e-03 CG   1
iter  8 act 9.416e-06 pre 7.823e-06 delta 4.186e-01 f 1.173e-05 |g| 8.019e-04 CG   1
[1] 0.961

结合pbdMPIRcpp Modules

让TRON呼叫分散式的运算系统

SPMD架构,无master,一次资料交换

objective_function <- function(w) {
  logger(sum(w))
  regularization <- sum((w - prior)^2) / 2
  lik <- sum(C * log(1 + exp(- y.value * Xv(X, w))))
  lik.all <- allreduce(x = lik, x.buffer = buffer.0, op = "sum")
  return(regularization + lik.all[1])
}

hs <- new(HsTrust, objective_function, ...)

SPMD架构

objective_function非常易于修改

所以我们能将注意力专注于模型上

objective_function <- function(w) {
  logger(sum(w))
  regularization <- sum((w - prior)^2) / 2
  lik <- sum(C * log(1 + exp(- y.value * Xv(X, w))))
  lik.all <- allreduce(x = lik, x.buffer = buffer.0, op = "sum")
  return(regularization + lik.all[1])
}

终于跨过资料量的门槛了...

该看看资料了!

资料越大,结果就会越好吗?

不同的模型对预测会有影响吗?

因子的组合

实验结果

auc Regularization FeatureSet DayEffect
1 1.0394 1 A 09:00:00
2 1.0233 10 A 09:00:00
3 1.0181 20 A 09:00:00
4 1.0326 1 A 16:00:00
5 1.0196 10 A 16:00:00
6 1.0125 20 A 16:00:00
7 1.0174 1 A 20:00:00
8 1.0057 10 A 20:00:00
9 1.0000 20 A 20:00:00
10 1.0436 1 B 09:00:00
11 1.0277 10 B 09:00:00
12 1.0215 20 B 09:00:00
13 1.0406 1 B 16:00:00
14 1.0269 10 B 16:00:00
15 1.0203 20 B 16:00:00
16 1.0217 1 B 20:00:00
17 1.0088 10 B 20:00:00
18 1.0017 20 B 20:00:00

分析

感谢R 强大的分析功能

row.names Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.0343 0.0009 1165.2725 0.0000
Regularization10 -0.0139 0.0009 -15.6506 0.0000
Regularization20 -0.0202 0.0009 -22.7578 0.0000
FeatureSetB 0.0049 0.0007 6.7933 0.0000
DayEffect20:00 -0.0162 0.0009 -18.2476 0.0000
DayEffect9:00 0.0035 0.0009 3.9816 0.0018

平均来说FeatureSet B 好 \(0.5\%\)

成果

plot of chunk unnamed-chunk-7

其他的工程面分享

处理大量数据需要很多工程

R Packages + Git + Jenkins:

开发后自动测试、部署和版本控制

利用Spark做Data Cleaning

R + rjson + awscli

自动利用AWS建立云端实验用Cluster

ec2_request_spot_instances <- function(spot_price, instance_count, launch_specification = list()) {
  write(toJSON(launch_specification), file=(json.path <- tempfile()))
  cmd <- sprintf("aws ec2 request-spot-instances --spot-price %f --instance-count %d --launch-specification %s", spot_price, instance_count, sprintf("file://%s", json.path))
  system_collapse(cmd)
}

R 是可以对大量数据进行处理:

使用Matrix套件的稀疏矩阵

Rcpp高效能的使用记忆体

Rcpp整合第三方的库

pbdMPI建立分散式平行运算丛集

关于今天分享的工作

感谢两位同学Y.-C. Juan, Y. Zhuang

和我在Bridgewell Inc.的合作

Q&A