转载

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

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

mxnet：结合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]