RでsvmをGPUで計算させてみた

とあることでビッグデータの解析でsvmを使うことになったのですが、あまりにも計算が遅くて使い物にならない・・・。で、サポートベクトルマシンの計算をGPUでできないか？と調べたところ、`Rgtsvm`というパッケージがあったので、試してみました。
前提となる環境は

- OS: Ubuntu 16.04 LTS
- R: Version 3.5.1
- CUDA: 9.0 (V9.0.176)
- NVIDIA Driver Version: 384.145
- libboost: 1.58

です

パッケージのソースは下記のGithubからダウンロードできます。
[Danko-Lab/Rgtsvm: The R package for SVM with GPU architecture based on the GTSVM software](https://github.com/Danko-Lab/Rgtsvm)
`e1071`パッケージとコンパチなのでパラメータの設定も同じ、グリッドサーチによるパラメータチューニングもできる、というところが非常によいです。使用する変数も`Matrix`,`SparseM`ライブラリによるスパース行列も扱えるみたいです。

#インストール
詳細はGithubに記載されているのでそちらを参照していただきたいのですが

- CUDAライブラリ
- Boostライブラリ
- `bit64`, `snow`, `SparseM`, `Matrix`のRパッケージ

が必要になります。GPUを使うのでCUDAのインストールは必須になります。インストール方法についてはここでは省略しますので、必要な方はググってください(苦笑

BoostライブラリについてはUbuntuであれば

```shell
$ sudo apt-get install libboost-all-dev
```

でインストールできますし、ソース[http://www.boost.org/users/download/#live](http://www.boost.org/users/download/#live)をダウンロードしてビルドして構いません。ビルド方法については[こちら](https://boostjp.github.io/howtobuild.html)が参考になります。

必要なライブラリが用意できたらビルドします。

```shell
$ git clone https://github.com/Danko-Lab/Rgtsvm.git
$ cd Rgtsvm
$ make R_dependencies
$ R CMD INSTALL --configure-args="--with-cuda-home=$YOUR_CUDA_HOME --with-boost-home=$YOUR_BOOST_HOME" Rgtsvm
```

`make R_dependencies`では必要なRライブラリ(`bit64`, `snow`, `SparseM`, `Matrix`)をインストールするのですが、既にインストール済みであれば省略してもOKです。repositoryの関係でダウンロードできない場合は、`rDeps.R`のrepoを変更or削除することで回避できますが、結局の所R上で`install.package`で個別にインストールしたほうが早いです。

ビルドの際は`--with-cuda-home=$YOUR_CUDA_HOME`にCUDAをインストールディレクトリを、`--with-boost-home=$YOUR_BOOST_HOME`にBoostをインストールディレクトリを指定します。

最初の方で述べたようにubuntu上で`apt-get`でlibboostをインストールした場合は、`/usr/local`を指定すると上手くいきます。

CUDA 9.0をインストールしている場合、ストリーミング マルチプロセッサ (SM) アーキテクチャがsm_20のGPUが使えないらしいです。GPUボードのアーキテクチャがどれなのかは[こちら](http://arnon.dk/matching-sm-architectures-arch-and-gencode-for-various-nvidia-cards/)で確認ができます。

デフォルトの設定ではアーキテクチャがsm_50になっていますが、GPUの環境に合わせて設定を変更することができます。例えばGTX-1080ならsm_60,Tesla V100ならsm_70という感じで。変更方法はビルド時の`--configure--args`に`--with-cuda-arch=sm_60`といった感じでアーキテクチャを指定します[^1]。変えることで何がどうなるか・・・はまでは検証していません。もしかするとデフォルトのままでもいいのかもしれません。

[^1]: 2021/2/10更新: 以前はconfigureファイルを直接変更していましたが仕様が変更されていたので記事を書き換えました。

<!--
デフォルトの設定ではアーキテクチャがsm_50になっていますが、GPUの環境に合わせて設定を変更することができます。例えばGTX-1080ならsm_60,Tesla V100ならsm_70という感じで。変更箇所は

```
# Rgtsvm\configure line 2381:
< NCFLAGS="-arch=sm_20 -O2"
---
> NCFLAGS="-arch=sm_60 -O2"

# Rgtsvm\configure.ac line 30:
< NCFLAGS="-arch=sm_20 -O2"
---
> NCFLAGS="-arch=sm_60 -O2"
```

です。変えることで何がどうなるか・・・はまでは検証していません。もしかするとデフォルトのままでもいいのかもしれません。
-->

#実行してみる
まずは`iris`で。

```r:R
require(Rgtsvm)
data(iris)
model <- svm(Species ~ ., data = iris)
print(model)

Call:
svm(formula = Species ~ ., data = iris)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  1 
      gamma:  0.25 
    tolerance:  0.001 
 time elapsed:  0.032 


Number of Support Vectors:  40
```

予測をしてみます。

```r:R
pred <- predict(model,newdata=iris)
table(pred, iris$Species)

pred         setosa versicolor virginica
  setosa         50          0         0
  versicolor      0         48         3
  virginica       0          2        47
```

と見慣れた結果画面が出てきます。が、あっという間に計算が終わってしまうのGPUが動いているかどうか確認できません。そこで、`Rgtsvm`の[vignette](https://github.com/Danko-Lab/Rgtsvm/blob/master/Rgtsvm-vignette.pdf)を参考にして、サンプルデータとして100000行×100列のビッグデータサンプルを作って計算させてみます。

<!--
(ちなみにこの[vignette](https://github.com/Danko-Lab/Rgtsvm/blob/master/Rgtsvm-vignette.pdf)ですが、あちこちにtypoがあるので注意が必要です)
-->

```r:R
set.seed(1)
size = 50000
dimension = 100
covar.mat <- matrix(runif(dimension * dimension), nrow = dimension)
covar.mat <- t(covar.mat) %*% covar.mat
zero <- mvrnorm(size, mu = c(1:dimension), Sigma = covar.mat)
one <- mvrnorm(size, mu = c(1:dimension) - 5, Sigma = covar.mat)
x <- rbind(zero, one)
y <- c(rep(0, nrow(zero)), rep(1, nrow(one)))
i.all <- 1:(2 * size)
i.training <- sample(i.all, length(i.all) * 0.8)
i.test <- i.all [!i.all %in% i.training]

tic()
model.gpu <-svm(x[i.training, ], y[i.training], type = "C-classification")
toc()

19.859 sec elasped
```
GPUにもよりますが、私が使用している環境(Tesla V100)ではおよそ20秒弱でした。

CPUではどうでしょうか？並列計算ができるsvmとして`parallelSVM`を使って計算してみます。使用するコアは20個です。

```r:R
library(parallelSVM)
tic()
model.cpu.parallel <- parallelSVM(
    x[i.training, ],
    y[i.training],
    type = "C-classification",
    numberCores = 20
    )
toc()

418.796 sec elasped
```
とおよそ20倍の速度に！シングルコアだと45分近くかかりました。約140倍!

予測計算も`Rgtsvm`はGPUを使いますので高速です。

```r:R
tic()
y.pred <-predict( model.gpu, x[i.test,] );
toc()

1.108 sec elasped

cat("accuracy", sum(y.pred == y[i.test])/length(i.test),"\n");

accuracy 0.8993
```

`parallelSVM`では`predict`をシングルコアで計算しているようなので、予測にとても時間がかかります。

```r:R
tic()
y.pred <-predict(model.cpu.parallel, x[i.test,] );
toc()

415.645 sec elasped

cat("accuracy", sum(y.pred == y[i.test])/length(i.test),"\n");

accuracy 0.89845 
```

なんと375倍!?

ちなみにGPUが動いているかどうかは、計算中に別窓でコンソールを開いて`nvidia-smi`を実行するとわかります。詳細な使い方はググってください。それほどGPU上のメモリを使っていないようです。

#パラメータについて
`Rgtsvm`は基本的に`e1071`ライブラリの`svm`と同じなのですが、いくつか`Rgtsvm`固有のパラメータがあります。

|Parameter　 |説明|
|:--------:|:----|
|type      |今の所は"C-classification"と"eps-regression"のみ|
|gpu.id    |マルチGPU環境の場合は計算に使用するGPUのID番号を指定。デフォルトは0|
|maxIter|最大繰り返し回数。defaultはNROW(x)*NCOL(x)|
|verbose|計算の途中経過を表示する。defaultはFALSE|

#パラメーターチューニング
`tune.svm`でパラメータチューニングも可能です。ただしマルチコアでは動かないので時間はかかります

```r:R
tic()
gt.tune <- tune.svm( x[i.training,], y[i.training],
    gamma = 2^seq(-11, -1, 2),
    cost = 10^(-1:1),
    tunecontrol=tune.control(sampling = "cross", cross=8, rough.cross=3)
    )
toc()

1237.995 sec elapsed

model<-svm( x[i.training,], y[i.training],
    gamma=as.numeric(gt.tune$best.parameters[1]),
    cost=as.numeric(gt.tune$best.parameters[2])
    )

y.pred <-predict.gtsvm(model, x[i.test,] );
cat("accuracy", sum(y.pred == y[i.test])/length(i.test),"\n");

accuracy 0.9252 
```

ということで、これまでは計算速度の関係でビッグデータでのSVM計算を遠慮していたのですが、`Rgtsvm`のおかげでいろいろと解析の幅が広がりそうです。