# Large Scale Learning in R

##

Wush Wu
Taiwan R User Group

Taiwan R User Group

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

主题包含但不限于：

R 套件使用

机器学习和统计模型

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

但是只要用对R 的包

R 是可以处理大量的数据

今天我想跟大家分享

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

主要是靠着`Rcpp`和`pbdMPI`等包

打造可扩放性的学习模组

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

## 实际的数据

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

## 课本的数据

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

今天从这样的数据开始

| 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`是pbdR(Programming with Big Data in R)专案的套件之一

## 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

### Lin et al. 2007 Trust Region Newton Methods for Large-scale Logistic Regression

为了加强效能

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

我们自己包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
- `grad`是`fun`的gradient
- `Hv`是`fun`的hessian乘上一个向量

## 利用`Rcpp Modules`把`function`暴露到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

结合`pbdMPI`和`Rcpp 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\%\)

成果

## 其他的工程面分享

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

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

- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50