所有由严酷的魔王发布的文章

mxnet:结合R与GPU加速深度学习

近年来,深度学习可谓是机器学习方向的明星概念,不同的模型分别在图像处理与自然语言处理等任务中取得了前所未有的好成绩。在实际的应用中,大家除了关心模型的准确度,还常常希望能比较快速地完成模型的训练。一个常用的加速手段便是将模型放在GPU上进行训练。然而由于种种原因,R语言似乎缺少一个能够在GPU上训练深度学习模型的程序包。

DMLC(Distributed (Deep) Machine Learning Community)是由一群极客发起的组织,主要目标是提供快速高质量的开源机器学习工具。近来流行的boosting模型xgboost便是出自这个组织。最近DMLC开源了一个深度学习工具mxnet,这个工具含有R,python,julia等语言的接口。本文以R接口为主,向大家介绍这个工具的性能与使用方法。

一、五分钟入门指南

在这一节里,我们在一个样例数据上介绍mxnet的基本使用方法。目前mxnet还没有登录CRAN的计划,所以安装方法要稍微复杂一些。

  • 如果你是Windows/Mac用户,那么可以通过下面的代码安装预编译的版本。这个版本会每周进行预编译,不过为了保证兼容性,只能使用CPU训练模型。
install.packages("drat", repos="https://cran.rstudio.com")
drat:::addRepo("dmlc")
install.packages("mxnet")
  • 如果你是Linux用户或者想尝试GPU版本,请参考这个链接里的详细编译教程在本地进行编译。

安装完毕之后,我们就可以开始训练模型了,下面两个小节分别介绍两种不同的训练神经网络的方法。

二分类模型与mx.mlp

首先,我们准备一份数据,并进行简单的预处理:

require(mlbench)
require(mxnet)
data(Sonar, package="mlbench")
Sonar[,61] = as.numeric(Sonar[,61])-1
train.ind = c(1:50, 100:150)
train.x = data.matrix(Sonar[train.ind, 1:60])
train.y = Sonar[train.ind, 61]
test.x = data.matrix(Sonar[-train.ind, 1:60])
test.y = Sonar[-train.ind, 61]

我们借用mlbench包中的一个二分类数据,并且将它分成训练集和测试集。mxnet提供了一个训练多层神经网络的函数mx.mlp,我们额可以通过它来训练一个神经网络模型。下面是mx.mlp中的部分参数:

  • 训练数据与预测变量
  • 每个隐藏层的大小
  • 输出层的结点数
  • 激活函数类型
  • 损失函数类型
  • 进行训练的硬件(CPU还是GPU)
  • 其他传给mx.model.FeedForward.create的高级参数

了解了大致参数后,我们就可以理解并让R运行下面的代码进行训练了。

mx.set.seed(0)
model <- mx.mlp(train.x, train.y, hidden_node=10, out_node=2,      out_activation="softmax", num.round=20, array.batch.size=15, learning.rate=0.07, momentum=0.9, eval.metric=mx.metric.accuracy)
## Auto detect layout of input matrix, use rowmajor..
## Start training with 1 devices
## [1] Train-accuracy=0.488888888888889
## [2] Train-accuracy=0.514285714285714
## [3] Train-accuracy=0.514285714285714

...

## [18] Train-accuracy=0.838095238095238
## [19] Train-accuracy=0.838095238095238
## [20] Train-accuracy=0.838095238095238

这里要注意使用mx.set.seed而不是R自带的set.seed函数来控制随机数。因为mxnet的训练过程可能会运行在不同的运算硬件上,我们需要一个足够快的随机数生成器来管理整个随机数生成的过程。模型训练好之后,我们可以很简单地进行预测:

preds = predict(model, test.x)
## Auto detect layout of input matrix, use rowmajor..
pred.label = max.col(t(preds))-1
table(pred.label, test.y)
##           test.y
## pred.label  0  1
##          0 24 14
##          1 36 33

如果进行的是多分类预测,mxnet的输出格式是类数X样本数。

回归模型与自定义神经网络

mx.mlp接口固然很方便,但是神经网络的一大特点便是它的灵活性,不同的结构可能有着完全不同的特性。mxnet的亮点之一便是它赋予了用户极大的自由度,从而可以任意定义需要的神经网络结构。我们在这一节用一个简单的回归任务介绍相关的语法。

首先,我们仍然要准备好一份数据。

data(BostonHousing, package="mlbench")

train.ind = seq(1, 506, 3)
train.x = data.matrix(BostonHousing[train.ind, -14])
train.y = BostonHousing[train.ind, 14]
test.x = data.matrix(BostonHousing[-train.ind, -14])
test.y = BostonHousing[-train.ind, 14]

mxnet提供了一个叫做“Symbol”的系统,从而使我们可以定义结点之间的连接方式与激活函数等参数。下面是一个定义没有隐藏层神经网络的简单例子:

# 定义输入数据
data <- mx.symbol.Variable("data")
# 完整连接的隐藏层
# data: 输入源
# num_hidden: 该层的节点数
fc1 <- mx.symbol.FullyConnected(data, num_hidden=1)

# 针对回归任务,定义损失函数
lro <- mx.symbol.LinearRegressionOutput(fc1)

在神经网络中,回归与分类的差别主要在于输出层的损失函数。这里我们使用了平方误差来训练模型。希望能更进一步了解Symbol的读者可以继续阅读这份以代码为主的文档。

定义了神经网络之后,我们便可以使用mx.model.FeedForward.create进行训练了。

mx.set.seed(0)
model <- mx.model.FeedForward.create(lro, X=train.x, y=train.y, ctx=mx.cpu(), num.round=50, array.batch.size=20, learning.rate=2e-6, momentum=0.9, eval.metric=mx.metric.rmse)
## Auto detect layout of input matrix, use rowmajor..
## Start training with 1 devices
## [1] Train-rmse=16.063282524034
## [2] Train-rmse=12.2792375712573
## [3] Train-rmse=11.1984634005885

...

## [48] Train-rmse=8.26890902770415
## [49] Train-rmse=8.25728089053853
## [50] Train-rmse=8.24580511500735

这里我们还针对回归任务修改了eval.metric参数。目前我们提供的评价函数包括”accuracy”,”rmse”,”mae” 和 “rmsle”,用户也可以针对需要自定义评价函数,例如:

demo.metric.mae <- mx.metric.custom("mae", function(label, pred) {
  res <- mean(abs(label-pred))
  return(res)
})
mx.set.seed(0)
model <- mx.model.FeedForward.create(lro, X=train.x, y=train.y, ctx=mx.cpu(), num.round=50, array.batch.size=20, learning.rate=2e-6, momentum=0.9, eval.metric=demo.metric.mae)
## Auto detect layout of input matrix, use rowmajor..
## Start training with 1 devices
## [1] Train-mae=13.1889538083225
## [2] Train-mae=9.81431959337658
## [3] Train-mae=9.21576419870059

...

## [48] Train-mae=6.41731406417158
## [49] Train-mae=6.41011292926139
## [50] Train-mae=6.40312503493494

至此,你已经掌握了基本的mxnet使用方法。接下来,我们将介绍更好玩的应用。

二、手写数字竞赛

在这一节里,我们以Kaggle上的手写数字数据集(MNIST)竞赛为例子,介绍如何通过mxnet定义一个强大的神经网络,并在GPU上快速训练模型。

第一步,我们从Kaggle上下载数据,并将它们放入data/文件夹中。然后我们读入数据,并做一些预处理工作。

require(mxnet)
train <- read.csv('data/train.csv', header=TRUE) 
test <- read.csv('data/test.csv', header=TRUE) 
train <- data.matrix(train) 
test <- data.matrix(test) 

train.x <- train[,-1] 
train.y <- train[,1]

train.x <- t(train.x/255)
test <- t(test/255)

最后两行预处理的作用有两个:

  • 原始灰度图片数值处在[0,255]之间,我们将其变换到[0,1]之间。
  • mxnet接受 像素X图片 的输入格式,所以我们对输入矩阵进行了转置。

接下来我们定义一个特别的神经网络结构:LeNet。这是Yann LeCun提出用于识别手写数字的结构,也是最早的卷积神经网络之一。同样的,我们使用Symbol语法来定义,不过这次结构会比较复杂。

# input
data <- mx.symbol.Variable('data')
# first conv
conv1 <- mx.symbol.Convolution(data=data, kernel=c(5,5), num_filter=20)
tanh1 <- mx.symbol.Activation(data=conv1, act_type="tanh")
pool1 <- mx.symbol.Pooling(data=tanh1, pool_type="max",
                          kernel=c(2,2), stride=c(2,2))
# second conv
conv2 <- mx.symbol.Convolution(data=pool1, kernel=c(5,5), num_filter=50)
tanh2 <- mx.symbol.Activation(data=conv2, act_type="tanh")
pool2 <- mx.symbol.Pooling(data=tanh2, pool_type="max",
                          kernel=c(2,2), stride=c(2,2))
# first fullc
flatten <- mx.symbol.Flatten(data=pool2)
fc1 <- mx.symbol.FullyConnected(data=flatten, num_hidden=500)
tanh3 <- mx.symbol.Activation(data=fc1, act_type="tanh")
# second fullc
fc2 <- mx.symbol.FullyConnected(data=tanh3, num_hidden=10)
# loss
lenet <- mx.symbol.SoftmaxOutput(data=fc2)

为了让输入数据的格式能对应LeNet,我们要将数据变成R中的array格式:

train.array <- train.x
dim(train.array) <- c(28, 28, 1, ncol(train.x))
test.array <- test
dim(test.array) <- c(28, 28, 1, ncol(test))

接下来我们将要分别使用CPU和GPU来训练这个模型,从而展现不同的训练效率。

n.gpu <- 1
device.cpu <- mx.cpu()
device.gpu <- lapply(0:(n.gpu-1), function(i) {
  mx.gpu(i)
})

我们可以将GPU的每个核以list的格式传递进去,如果有BLAS等自带矩阵运算并行的库存在,则没必要对CPU这么做了。

我们先在CPU上进行训练,这次我们只进行一次迭代:

mx.set.seed(0)
tic <- proc.time()
model <- mx.model.FeedForward.create(lenet, X=train.array, y=train.y, ctx=device.cpu, num.round=1, array.batch.size=100, learning.rate=0.05, momentum=0.9, wd=0.00001, eval.metric=mx.metric.accuracy, epoch.end.callback=mx.callback.log.train.metric(100))
## Start training with 1 devices
## Batch [100] Train-accuracy=0.1066
## Batch [200] Train-accuracy=0.16495
## Batch [300] Train-accuracy=0.401766666666667
## Batch [400] Train-accuracy=0.537675
## [1] Train-accuracy=0.557136038186157
print(proc.time() - tic)
##    user  system elapsed
## 130.030 204.976  83.821

在CPU上训练一次迭代一共花了83秒。接下来我们在GPU上训练5次迭代:

mx.set.seed(0)
tic <- proc.time()
model <- mx.model.FeedForward.create(lenet, X=train.array, y=train.y, ctx=device.gpu, num.round=5, array.batch.size=100, learning.rate=0.05, momentum=0.9, wd=0.00001, eval.metric=mx.metric.accuracy, epoch.end.callback=mx.callback.log.train.metric(100))
## Start training with 1 devices
## Batch [100] Train-accuracy=0.1066
## Batch [200] Train-accuracy=0.1596
## Batch [300] Train-accuracy=0.3983
## Batch [400] Train-accuracy=0.533975
## [1] Train-accuracy=0.553532219570405
## Batch [100] Train-accuracy=0.958
## Batch [200] Train-accuracy=0.96155
## Batch [300] Train-accuracy=0.966100000000001
## Batch [400] Train-accuracy=0.968550000000003
## [2] Train-accuracy=0.969071428571432
## Batch [100] Train-accuracy=0.977
## Batch [200] Train-accuracy=0.97715
## Batch [300] Train-accuracy=0.979566666666668
## Batch [400] Train-accuracy=0.980900000000003
## [3] Train-accuracy=0.981309523809527
## Batch [100] Train-accuracy=0.9853
## Batch [200] Train-accuracy=0.985899999999999
## Batch [300] Train-accuracy=0.986966666666668
## Batch [400] Train-accuracy=0.988150000000002
## [4] Train-accuracy=0.988452380952384
## Batch [100] Train-accuracy=0.990199999999999
## Batch [200] Train-accuracy=0.98995
## Batch [300] Train-accuracy=0.990600000000001
## Batch [400] Train-accuracy=0.991325000000002
## [5] Train-accuracy=0.991523809523812
print(proc.time() - tic)
##    user  system elapsed
##   9.288   1.680   6.889

在GPU上训练5轮迭代只花了不到7秒,快了数十倍!可以看出,对于这样的网络结构,GPU的加速效果是非常显著的。有了快速训练的办法,我们便可以很快的做预测,并且提交到Kaggle上了:

preds <- predict(model, test.array)
pred.label <- max.col(t(preds)) - 1
submission <- data.frame(ImageId=1:ncol(test), Label=pred.label)
write.csv(submission, file='submission.csv', row.names=FALSE, quote=FALSE)

三、图像识别应用

其实对于神经网络当前的应用场景而言,识别手写数字已经不够看了。早些时候,Google公开了一个云API,让用户能够检测一幅图像里面的内容。现在我们提供一个教程,让大家能够自制一个图像识别的在线应用。

DMLC用在ImageNet数据集上训练了一个模型,能够直接拿来对真实图片进行分类。同时,我们搭建了一个Shiny应用,只需要不超过150行R代码就能够自己在浏览器中进行图像中的物体识别。

为了搭建这个应用,我们要安装shiny和imager两个R包:

install.packages("shiny", repos="https://cran.rstudio.com")
install.packages("imager", repos="https://cran.rstudio.com")

现在你已经配置好了mxnet, shiny和imager三个R包,最困难的部分已经完成了!下一步则是让shiny直接下载并运行我们准备好的代码:

shiny::runGitHub("thirdwing/mxnet_shiny")

第一次运行这个命令会花上几分钟时间下载预先训练好的模型。训练的模型是Inception-BatchNorm Network,如果读者对它感兴趣,可以阅读这篇文章。准备就绪之后,你的浏览器中会出现一个网页应用,就用本地或在线图片来挑战它吧!

如果你只需要一个图像识别的模块,那么我们下面给出最简单的一段R代码让你能进行图像识别。首先,我们要导入预训练过的模型文件:

model <<- mx.model.load("Inception/Inception_BN", iteration = 39)
synsets <<- readLines("Inception/synset.txt")
mean.img <<- as.array(mx.nd.load("Inception/mean_224.nd")[["mean_img"]])

接下来我们使用一个函数对图像进行预处理,这个步骤对于神经网络模型而言至关重要。

preproc.image <- function(im, mean.image) {
  # crop the image
  shape <- dim(im)
  short.edge <- min(shape[1:2])
  yy <- floor((shape[1] - short.edge) / 2) + 1
  yend <- yy + short.edge - 1
  xx <- floor((shape[2] - short.edge) / 2) + 1
  xend <- xx + short.edge - 1
  croped <- im[yy:yend, xx:xend,,]
  # resize to 224 x 224, needed by input of the model.
  resized <- resize(croped, 224, 224)
  # convert to array (x, y, channel)
  arr <- as.array(resized)
  dim(arr) = c(224, 224, 3)
  # substract the mean
  normed <- arr - mean.img
  # Reshape to format needed by mxnet (width, height, channel, num)
  dim(normed) <- c(224, 224, 3, 1)
  return(normed)
}

最后我们读入图像,预处理与预测就可以了。

im <- load.image(src)
normed <- preproc.image(im, mean.img)
prob <- predict(model, X = normed)
max.idx <- order(prob[,1], decreasing = TRUE)[1:5]
result <- synsets[max.idx]

四、参考资料

MXNet是一个在底层与接口都有着丰富功能的软件,如果读者对它感兴趣,可以参考一些额外的材料来进一步了解MXNet,或者是深度学习这个领域。

 

xgboost: 速度快效果好的boosting模型

本文作者:何通,SupStat Inc(总部在纽约,中国分部为北京数博思达信息科技有限公司)数据科学家,加拿大Simon Fraser University计算机学院研究生,研究兴趣为数据挖掘和生物信息学。

主页:https://github.com/hetong007

引言

在数据分析的过程中,我们经常需要对数据建模并做预测。在众多的选择中,randomForest, gbmglmnet是三个尤其流行的R包,它们在Kaggle的各大数据挖掘竞赛中的出现频率独占鳌头,被坊间人称为R数据挖掘包中的三驾马车。根据我的个人经验,gbm包比同样是使用树模型的randomForest包占用的内存更少,同时训练速度较快,尤其受到大家的喜爱。在python的机器学习库sklearn里也有GradientBoostingClassifier的存在。

Boosting分类器属于集成学习模型,它基本思想是把成百上千个分类准确率较低的树模型组合起来,成为一个准确率很高的模型。这个模型会不断地迭代,每次迭代就生成一颗新的树。对于如何在每一步生成合理的树,大家提出了很多的方法,我们这里简要介绍由Friedman提出的Gradient Boosting Machine。它在生成每一棵树的时候采用梯度下降的思想,以之前生成的所有树为基础,向着最小化给定目标函数的方向多走一步。在合理的参数设置下,我们往往要生成一定数量的树才能达到令人满意的准确率。在数据集较大较复杂的时候,我们可能需要几千次迭代运算,如果生成一个树模型需要几秒钟,那么这么多迭代的运算耗时,应该能让你专心地想静静……

Rocket

现在,我们希望能通过xgboost工具更好地解决这个问题。xgboost的全称是eXtreme Gradient Boosting。正如其名,它是Gradient Boosting Machine的一个c++实现,作者为正在华盛顿大学研究机器学习的大牛陈天奇。他在研究中深感自己受制于现有库的计算速度和精度,因此在一年前开始着手搭建xgboost项目,并在去年夏天逐渐成型。xgboost最大的特点在于,它能够自动利用CPU的多线程进行并行,同时在算法上加以改进提高了精度。它的处女秀是Kaggle的希格斯子信号识别竞赛,因为出众的效率与较高的预测准确度在比赛论坛中引起了参赛选手的广泛关注,在1700多支队伍的激烈竞争中占有一席之地。随着它在Kaggle社区知名度的提高,最近也有队伍借助xgboost在比赛中夺得第一

为了方便大家使用,陈天奇将xgboost封装成了python库。我有幸和他合作,制作了xgboost工具的R语言接口,并将其提交到了CRAN上。也有用户将其封装成了julia库。python和R接口的功能一直在不断更新,大家可以通过下文了解大致的功能,然后选择自己最熟悉的语言进行学习。

继续阅读xgboost: 速度快效果好的boosting模型

自制简单遗传算法实验

我参加完八月份的COS沙龙之后比较闲,忽然想起自己很久以前看的遗传算法的基本思想。本着时不时就应该做一些私活的心态,我就在旅行商问题上面把它实现了一下。

遗传算法

遗传算法是一个仿生学的算法。进化论认为地球上千奇百怪的生物都是进化而来的,如今能生存在地球上的生物是更适应于这个环境的,我们也可以说它们是被“优化过”的。他们是怎么优化的呢?在一个种群中,生物的差异主要来自于两点,不同染色体之间的交叉结合以及染色体自发的随机变异。这些差异实际上是随机发生的,但是生物的外部生存环境会通过生存与死亡让更能适应的个体存活下去。因此随着时间的推移,生物种群对环境的适应能力会越来越高。受到这个现象的启发,有人发明了遗传算法,通过模拟遗传的过程来解决一些优化问题。
继续阅读自制简单遗传算法实验

在R中实现动态气泡图

最近我逐渐发现了ggplot2这个包的好处——只要用过一次,就再也不想回头使用R中自带的作图函数了。前两天鼓捣完一个地图的数据,又受到统计之都最新文章的影响,我忽然想起了Hans Rosling在TED上的精彩演讲。在图中横坐标是国民收入,纵坐标是国民的期望寿命,气泡的大小则是该国人口。整个图从1800年的统计数据开始,一直到2009年不断动态地展示,图上的气泡也随着时间变化不停地抖动上升。有一位在斯坦福专做可视化的博士用JavaScript在网页上重现了这段动态效果图,点开页面即可观看:http://bost.ocks.org/mike/nations/

于是今天我便将这个图尝试着用R中的ggplot2与animation包实现了出来,边实现边研究ggplot2的用法,花了一天的时间做成了下面的这个视频

简单地说一下流程:首先是数据文件的获取。数据能够在github上找到,但是数据是JSON格式的,只有一行,因此我的大部分代码都在为让数据变成一个二维矩阵的形式而努力着……很多国家会出现某些年没有统计数据的情况,因此我用了线性插值填补。最后,有两个国家只有一年有数据,我只能将它们删掉了。

弄好了数据就可以使用ggplot2画图了。为了让图像好看,我调整了图像的属性,比如圆圈的大小范围,学习加边框,学习图中加文字(annotate)等语法。但我现在感觉还是有一些地方能够微调改进。

最后使用animation包中的saveMovie函数,结合ffmpeg导出成了一个视频。

最后附上代码