---
title       : Large Scale Learning in R
subtitle    : 
author      : Wush Wu
job         : Taiwan R User Group
framework   : io2012        # {io2012, html5slides, shower, dzslides, ...}
highlighter : highlight.js  # {highlight.js, prettify, highlight}
hitheme     : tomorrow       # 
widgets     : [mathjax]            # {mathjax, quiz, bootstrap}
mode        : selfcontained # {standalone, draft}
license: by-nc-sa
logo: Taiwan-R-logo.png
--- &vcenternobg .large

```{r setup, include=FALSE,echo = F, message = F, warning = F, tidy = F, cache=FALSE}
# make this an external chunk that can be included in any file
library(xtable)
library(rbenchmark)
library(reshape2)
library(knitr)
library(methods)
options(width = 100,digits=3)
opts_chunk$set(message = FALSE, eval=TRUE,fig.align = "center", warning = FALSE, comment = NA, dpi = 100, fig.width=6, fig.height=4.5,tidy = FALSE, cache = FALSE, echo=FALSE)

options(xtable.type = 'html')
knit_hooks$set(inline = function(x) {
  if(is.numeric(x)) {
    round(x, getOption('digits'))
  } else {
    paste(as.character(x), collapse = ', ')
  }
})
knit_hooks$set(plot = knitr:::hook_plot_html)
library(bibtex)
library(knitcitations)
bib <- read.bib("index.bib")
```

<img src="assets/img/MLDMMonday.jpg" class="grayscale fit100" />

Taiwan R User Group

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

主题包含但不限于：

R 套件使用

机器学习和统计模型

--- &vcenter .large

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

--- &vcenter .large

但是只要用对R 的包

R 是可以处理大量的数据

--- &vcenter .large

今天我想跟大家分享

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

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

打造可扩放性的学习模组

--- &vcenter .large

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

--- &twocolvcenter

*** =left

## 实际的数据

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

*** =right

## 课本的数据

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

--- &vcenter .large

今天从这样的数据开始

```{r demo_data, results='asis'}
df <- read.table("demo_data.txt")
colnames(df) <- c("is_click", "show_time", "", "client_ip")
df$adid <- rownames(df)
rownames(df) <- NULL
df$show_time <- paste(df[,2], df[,3])
df[[3]] <- NULL
print(xtable(df), type="html")
```
<br/>
以“预测点击发生的机率”为例

--- &twocol .large

## 问题的建模

```{r demo_data2, dependson="demo_data", results='asis'}
print(xtable(df), type="html")
```
<h2>$$\Downarrow$$</h2>
<h2>$$P(y|x) = \frac{1}{1 + e^{-yw^Tx}}$$</h2>

--- &vcenter .large

$$w^Tx$$

--- &vcenter .large

iris

```{r iris_way, echo=TRUE}
head(model.matrix(Species ~ ., iris))
```

--- &vcenter .large

实际的资料

有 $10^{9+}$ 笔资料

有 $10^{3+}$ 笔广告

每种类别变量可能到 $10^{3}$ 类别

--- &vcenter .large

`model.matrix`的后果：

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

也就是 $`r 10^12 / 2^30 / 8`$ GB的记忆体

--- &vcenter .large

开始抱怨了!

### R 不能处理大量的数据

### R 会吃很多记忆体

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

--- &vcenter .large

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

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

--- &twocolvcenter

*** =left

## 稀疏矩阵

- 只储存非0的资讯

```{r sparse_matrix_echo, echo=TRUE, eval=FALSE}
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)
```

*** =right

```{r sparse_matrix_eval, cache=TRUE}
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)
```

--- &vcenter .large

在我手上的资料上

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

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

--- &vcenter .large

如果`m1`, `m2`是 $10^4 \times 10^4$ 的矩阵:

```{r sparse_mat_efficiency, cache=TRUE}
library(rbenchmark)
n <- 4
m1 <- matrix(0, 10^n, 10^n);m1[1, 4] <- 1
m2 <- Matrix(0, 10^n, 10^n, sparse=TRUE)
r <- rnorm(10^n)
benchmark(
  m1 %*% r,
  m2 %*% r)
```

--- &vcenter .large

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

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

--- &vcenter .large

```{r rcpp_matrix_multiplication_eval, cache=TRUE, engine='Rcpp'}
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
SEXP XTv(S4 m, NumericVector v, NumericVector& retval) {
  int
    *i = INTEGER(wrap(m.slot("i"))),
    i_len = Rf_length(wrap(m.slot("i"))),
    *p = INTEGER(wrap(m.slot("p"))),
    p_len = Rf_length(wrap(m.slot("p")));
    NumericVector x = wrap(m.slot("x"));
  IntegerVector dim(wrap(m.slot("Dim")));
  if (v.size() != dim[0]) throw std::invalid_argument("");
  for(int j = 0,k = 0;j < p_len - 1;j++) {
    while(k < p[j+1]) {
      retval[j] += x[k]*v[i[k++]];
    }
  }
  return retval;
}
```

```{r rcpp_matrix_multiplication_benchmark, cache=TRUE, dependson=c("sparse_mat_efficiency","rcpp_matrix_multiplication_eval"), echo=FALSE}
retval <- numeric(10^n)
benchmark(
  m2 %*% r,
  XTv(m2, r, retval))
```

在现代的Rcpp架构下

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

我们只需要专注在算法上

```{r rcpp_matrix_multiplication_echo, cache=TRUE, engine='Rcpp', echo=TRUE}
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
SEXP XTv(S4 m, NumericVector v, NumericVector& retval) {
  //...
}
```

--- &vcenter .large

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

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

--- &vcenter .large

$$\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$$

--- &twocolvcenter .large

***=left

## MPI

### 记忆体足够的话，较快

### 没有容错

***=right

## Hadoop

### 慢，要大量机器才有效果

### 有容错

--- &vcenter .large

`pbdMPI`是[pbdR(Programming with Big Data in R)](http://r-pbd.org/)专案的套件之一

<img src="assets/img/pbdMPI.png" class="fit50" />

--- &twocolvcenter .large

*** =left

## pbdMPI的优点

### 好安装
### 好开发

*** =right

```{r pbdMPI_demo, eval=FALSE, echo=TRUE, tidy.opts=list(width.cutoff=20L)}
library(pbdMPI)
.rank <- comm.rank()
filename <- sprintf("%d.csv", .rank)
data <- read.csv(filename)
target <- reduce(sum(data$value), op="sum")
finalize()
```

--- &twocolvcenter .large

*** =left

## 最佳化算法

### 迭代次数少
### 不用计算Hessian矩阵
### 已有很棒的实作[LIBLINEAR](http://www.csie.ntu.edu.tw/~cjlin/liblinear/)

*** =right

<h3>`r citet(bib["Lin:2007:TRN:1273496.1273567"], format_inline_fn = function(x) paste(knitcitations:::format_authoryear_p(x), x$title))`</h3>

--- &vcenter .large

为了加强效能

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

我们自己包LIBLINEAR到R中

--- &vcenter .large

## LIBLINEAR 原始码中的 `tron.h`

```{r tron, engine='Rcpp', eval=FALSE, echo=TRUE}
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`提供了一个界面

--- &vcenter

## 界面：`function`

```{r function_interface, engine='Rcpp', eval=FALSE, echo=TRUE}
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乘上一个向量

--- &twocolvcenter .large

## 利用`Rcpp Modules`把`function`暴露到R中

*** =left

### 实作Rfunction

```{r Rfunction, engine='Rcpp', eval=FALSE, echo=TRUE}
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//...
```

*** =right

将Rfunction暴露出来

```{r Rfunction_module, engine='Rcpp', eval=FALSE, echo=TRUE}
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)
	;	
}
```

--- &twocolvcenter .large

*** =left

实作的结果可包装成[套件](https://bitbucket.org/wush978/largescalelogisticregression/src/4daf9c5bba5cd0e4f35afd813866e6da72ca92bb/?at=hstrust)


```{r, eval=FALSE, echo=TRUE}
library(devtools)
install_bitbucket(
  repo="largescalelogisticregression", 
  username="wush978", ref="hstrust")
```

*** =right

使用

```{r, eval = FALSE, echo=TRUE}
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))
```

--- &vcenter

```{r, eval = FALSE, echo=TRUE}
library(HsTrust)
# ...
print(r <- obj$tron(1e-4, TRUE))
```

```{r, eval = TRUE, echo=FALSE}
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))
```


--- &vcenter .large

结合`pbdMPI`和`Rcpp Modules`

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

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

```{r, eval=FALSE, echo=TRUE}
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, ...)
```

--- &vcenter .large

SPMD架构

<img src="assets/img/SPMD.png" class="fit100"/>


--- &vcenter .large

`objective_function`非常易于修改

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

```{r, eval=FALSE, echo=TRUE}
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])
}
```

--- &vcenter .large

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

该看看资料了！

--- .dark .segue

## 资料越大，结果就会越好吗？

---

<img src="assets/img/datasize.png" class="fit100"/>

--- .dark .segue

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

--- &vcenter .large

因子的组合

<img src="assets/img/factors.png" class="fit100" />

--- &vcenter .large

实验结果

```{r exp, echo=FALSE, results='asis'}
print(xtable(head(read.csv("exp.csv", header=TRUE), 20), digits=4), type="html")
```

--- &vcenter .large

分析

### 感谢R 强大的分析功能

```{r result, echo=FALSE, results='asis'}
print(xtable(read.csv("result.csv", sep="\t", check.names=FALSE), digits=4), type="html")
```

### 平均来说FeatureSet B 好 $0.5\%$

--- &vcenter .large

<img src="assets/img/cycle.png" class="fit100" />

--- &vcenter .large

成果

```{r, echo=FALSE}
library(ggplot2)
df2 <- readRDS("../result.Rds")
ggplot(df2, aes(x=date, y=ratio)) +
  geom_line()
```

--- .dark .segue

## 其他的工程面分享

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

--- &vcenter .large

R Packages + Git + Jenkins:

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

<img src="assets/img/jenkins.png" class="fit50"/>

--- &vcenter .large

利用Spark做Data Cleaning

<img src="assets/img/Spark.png" class="fit50"/>

--- &vcenter .large

R + [rjson]() + [awscli](https://aws.amazon.com/cn/cli/)

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

```{r, eval=FALSE, echo=TRUE}
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)
}
```

--- &vcenter .large

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

### 使用`Matrix`套件的稀疏矩阵
### `Rcpp`高效能的使用记忆体
### `Rcpp`整合第三方的库
### `pbdMPI`建立分散式平行运算丛集

--- &vcenter .large

关于今天分享的工作

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

和我在Bridgewell Inc.的合作

--- &vcenter .large

Q&A
