はじめに
どうも，Machine Learning Advent Calendar 2012の3日目を担当させていただきます．Advent Calendarに参加させていただくのは初めてなので恐縮ですが，よろしくお願いいたします．本日の内容は，githubで公開しました．
本日の内容
さて，本日の内容ですが，タイトルの通り，OpenMPを使ってLIBSVMをマルチコア環境で速くする方法です．LIBSVMは，多分皆さんご存知だと思いますが，Support Vector Machines (SVM)の代表的で信頼出来る実装です．台湾大学のChih-Chung Chang，Chih-Jen Linらによって開発されています．（直接関係ありませんが，台湾は中国語のローマ字表記に，有名なピンインではなくウェード式を用いているため，"Chih"のような普通のピンイン表記では現れない綴りが名前の表記に使われています．）
SVMの特徴の１つは，カーネルトリックを使って，複雑な識別面を構成できることです．まぁ，複雑な数式はおいておくとして，LIBSVMではオプションを付けるだけなのですが，実用上は，カーネルを使うと（特にデータが多い場合は）学習が遅い，という問題があります．
…学習，速くしたいですよね？（笑）せっかく，コアが沢山あるのですから．実は，公式のFAQでOpenMPを使って並列化する方法が公開されているのです！早速，このLIBSVM FAQを見てましょう．
基本的には，Makefileに-fopenmpを付け加えて，svm.cppに3行付け加えるだけです．一応，現時点での最新バージョン3.14用のパッチファイルを作ってみましたが，パッチファイルはLIBSVMがバージョンアップしたら使えません．このFAQ，ちゃんとアップデートされていますので，しばらく経ってからこのページを見る方は，直接，LIBSVM FAQを参照して下さい．
wget 'http://www.csie.ntu.edu.tw/~cjlin/cgi-bin/libsvm.cgi?+http://www.csie.ntu.edu.tw/~cjlin/libsvm+tar.gz'
git clone https://github.com/niam/libsvm_openmp_patch
tar zxvf libsvm-3.14.tar.gz
cd libsvm-3.14
patch < ../libsvm_openmp_patch/Makefile.diff
patch < ../libsvm_openmp_patch/svm_cpp.diff
make
コンパイルには，OpenMPが動作する環境が必要です．確認は，gcc 4.7.2で行なっています．OpenMPについ簡単に説明すると，for文の前に#pragma omp parallel forと書くと，そのfor文を並列にしてくれる，超お手軽並列化です．実際に何並列にするかは，OMP_NUM_THREADS環境変数で動的に設定できます．どれぐらい速くなるのか，実験してみました．テストのやり方も，先程のFAQに書いてあるので，それに従います．
wget 'http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary/ijcnn1.bz2'
bzip2 -d ijcnn1.bz2
まずは，1スレッドの時．
export OMP_NUM_THREADS=1
time svm-train -c 16 -g 4 -m 400 -t 2 ijcnn1
とします．OMP_NUM_THREADSで，スレッド数を設定しています．-m 400は，カーネルの値をキャッシュしておくサイズをMB単位で指定するもので，ここでは400MB使っています．-t 2は，RBFカーネルを使用する事を表し，-g 4はRBFカーネルのパラメータです．-c 16は正則化の度合いです．精度の評価をする時などは，-g, -cなどの値を，交差検定などを用いて決めてやらなければなりませんが，ここでは速度だけを見たいので，決め打ちにしています．
...............*.........*.*
optimization finished, #iter = 25062
nu = 0.025946
obj = -14014.950327, rho = -0.415715
nSV = 3370, nBSV = 792
Total nSV = 3370
./svm-train -c 16 -g 4 -t 2 -m 400 ijcnn1 38.67s user 0.21s system 99% cpu 38.889 total
でした．一方，8スレッドの時は…
export OMP_NUM_THREADS=8
time svm-train -c 16 -g 4 -t 2 -m 400 ijcnn1
とすると，
...............*.........*.*
optimization finished, #iter = 25062
nu = 0.025946
obj = -14014.950327, rho = -0.415715
nSV = 3370, nBSV = 792
Total nSV = 3370
./svm-train -c 16 -g 4 -t 2 -m 400 ijcnn1 55.77s user 0.35s system 559% cpu 10.038 total
38.889秒→10.038秒なので，約4倍速くなってますね！
僕の研究室のマシンの１つは，24コアあります．これをフルに使うと…？
export OMP_NUM_THREADS=24
time svm-train -c 16 -g 4 -t 2 -m 400 ijcnn1
とすると，
...............*.........*.*
optimization finished, #iter = 25062
nu = 0.025946
obj = -14014.950327, rho = -0.415715
nSV = 3370, nBSV = 792
Total nSV = 3370
./svm-train -c 16 -g 4 -t 2 -m 400 ijcnn1 131.44s user 0.97s system 1463% cpu 9.047 total
でした！
おまけ
LIBSVMは，様々なカーネル関数を選択するのに関数ポインタを使った動的多態を使っています．カーネル関数は，データの数をnとすると，O(n^2)回も呼び出される関数です．C++にちょっと詳しい人は，「え，それ，呼び出しコスト高くね？静的多態を使ってみた方が良くね？」と思うかも知れません．なので，呼び出しコストがどれぐらいか，ちょっと試してみました．
#include <iostream>
#include <functional>
#include <memory>
#include <chrono>
#include <cmath>
struct rbf_t{
inline double operator()(int i, int j) const{
double tmp = i-j;
return 1.0;
// return std::exp(-tmp*tmp);
}
};
struct kbase_t{
virtual double operator()(int i, int j) const =0;
};
struct rbf_t2:kbase_t{
inline double operator()(int i, int j) const{
double tmp = i-j;
return 1.0;
// return std::exp(-tmp*tmp);
}
};
inline double dyn_rbf(int i, int j) {
double tmp = i-j;
return 1.0;
// return std::exp(-tmp*tmp);
}
template<class K>
void kernelsum(K k, int n, const std::string& str){
auto time_point_t1 = std::chrono::high_resolution_clock::now();
double sum = 0.;
for(int i=0;i<n;++i){
for(int j=0;j<n;++j){
sum += k(i,j);
}
}
auto duration1 = std::chrono::high_resolution_clock::now() - time_point_t1 ;
std::cout << str << " time: " << duration1.count() << std::endl;
std::cout << str << " val : " << sum << std::endl;
return;
}
void kernelsum2(const kbase_t& k, int n, const std::string& str){
auto time_point_t1 = std::chrono::high_resolution_clock::now();
double sum = 0.;
for(int i=0;i<n;++i){
for(int j=0;j<n;++j){
sum += k(i,j);
}
}
auto duration1 = std::chrono::high_resolution_clock::now() - time_point_t1 ;
std::cout << str << " time: " << duration1.count() << std::endl;
std::cout << str << " val : " << sum << std::endl;
return;
}
int main(int argc, char* argv[]){
int N = 100000;
kernelsum(rbf_t(), N, std::string("static polymorphism "));
rbf_t2 rbf2;
kernelsum2(rbf2, N, std::string("dynamic polymorphism "));
std::function<double ( int, int)> f = &dyn_rbf;
kernelsum(f, N, std::string("dynamic polymorphism using std::function "));
}
関数の呼び出しコストを計測するだけなので，カーネル関数は，全ての場合に1.0を返す簡単なものになっています．C++11では，関数や関数オブジェクトをまとめてstd::functionに突っ込むことが出来ます．また，std::chronoを使って，実行時間ををμs単位で計測することができます．コンパイルは：
g++ -O3 --std=c++11 -o mykernel_test2 mykernel_test2.cpp
実行すると：
static polymorphism time: 8749249
static polymorphism val : 1e+10
dynamic polymorphism time: 29043178
dynamic polymorphism val : 1e+10
dynamic polymorphism using std::function time: 26236452
dynamic polymorphism using std::function val : 1e+10
単位はμsです．とりあえず，データ数が100,000の場合を想定すると，静的多態が8.7秒ぐらい関数呼び出しにかかっているのに対して，仮想関数を使った動的多態は29秒，std::functionを使った動的多態は26秒かかっていることがわかります．データ数が100,000の場合に，現実的なデータに対してRBFカーネルのSVMを学習しようとすると，普通は数時間以上はかかるはずです．なので，簡単にいえば，学習にかかる時間の方がずっと多いため，関数呼び出しコストをそんなに気にする必要はない，という結論になりました．
終わりに
本当は，SVMの式を出して，カーネル関数の呼び出しの部分を並列化しているところまでちゃんと解説したかったのですが，githubにpushできないエラーにハマって（よくわからないですが解決しました）心が折れてしまったので，この辺で終わりにしたいと思います．
最後に，主催者で，このAdvent CalendarのためにQiitaでTeXを扱うextensionまで作って下さったnaoya_tさんに感謝したいと思います．Chromeで，拡張機能設定画面に.crxファイルをドラッグ＆ドロップすると，extensionがインストールできる事を初めて知りました…．