C++11 で線形分類器の分散オンライン学習器を実装してみた

自作の（非）線形分類器のオンライン学習器に，

• Distributed Training Strategies for the Structured Perceptron (NAACL 2010)

で提案されている分散オンライン学習法 (Iterative Parameter Mixing) を実装してみた．

• 機械学習 × MapReduce - ny23の日記

で調べたのがもう2年以上前のことだから今更感はんぱない．と言っても，自作の分類器に実装したのは線形学習の分散化で，多項式カーネルを使った非線形学習の分散化はパラメタを効率的に分配するのが難しく，まだ実装できていない．というか，線形学習の分散化も平均化を行う場合はパラメタの管理が面倒で，効率的な実装になかなか辿りつけなかった．
マルチコア CPU 向けの分散オンライン学習は C++11 で新たに導入された thread と lambda 式を使うことで，驚くほど簡単に実装できる．既存の実装 +10~15 行程度で済む，ということを示すために書いたコード (分散 Passive Aggressive-I アルゴリズム) は以下．該当箇所は，distributed training と書かれた行から non-distributed training までで，分散化に必要となるパラメタを合わせると14行．

// ppa1.cc ; GNU GPL version 2 copyright@id:ny23
#include <err.h>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <chrono>
#include <thread>

static const size_t MAX_LINE_LEN = 1 << 20; // max line len
static float power10[40];

typedef std::vector <std::pair <int, float> >  x_t;
typedef std::pair <int, x_t>  ex_t;
struct result_t { int pp, pn, np, nn; };

int strtoi (char*& p) {
#ifdef USE_MY_STRTOI // faster strtoi
  int i = 0;
  for (; *p >= '0' && *p <= '9'; ++p) i *= 10, i += *p, i -= '0';
  return i;
#else
  return static_cast <int> (std::strtol (p, &p, 10));
#endif
}

float strtof (char*& p) {
#ifdef USE_MY_STRTOF // faster strtof; only support floating numbers -?\d*.\d*
  float i = 0; float sign = 1.0f; char* q = 0;
  if (*p == '-') ++p, sign = -1.0f;
  while ((*p >= '0' && *p <= '9') || *p == '.')
    if (*p == '.') q = ++p; else i *= 10, i += *p, i -= '0', ++p;
  return q ? sign * i * power10[p - q] : sign * i;
#else
  return static_cast <float> (std::strtod (p, &p)); // std::strtof is slower
#endif
}

// read examples
void read_data (const char* train, std::vector <ex_t>& ex, int* imax) {
  FILE* reader = std::fopen (train, "r");
  char line[MAX_LINE_LEN];
  for (x_t x; std::fgets (line, MAX_LINE_LEN, reader) != NULL; x.clear ()) {
    char* p = line;
    const int y = static_cast <int> (std::strtol (p, &p, 10));
    float norm = 0;
    while (*++p) {
      const int   i = strtoi (p);
      const float v = strtof (++p);
      norm += v * v;
      x.emplace_back (i, v);
    }
    norm = std::sqrt (norm);
    for (auto& f : x) f.second /= norm;
    if (! x.empty ()) *imax = std::max (*imax, x.back ().first);
    ex.emplace_back (y, x);
  }
  std::fclose (reader);
}

// one epoch passive aggressive I
void train1 (const ex_t* const start, const ex_t* const end, std::vector <float>& w, const float c) {
  for (const ex_t* p = start; p != end; ++p) {
    float m = 0;
    for (const auto& f : p->second) m += w[f.first] * f.second;
    m *= p->first;
    if (m <= 1.0f) {
      const float t = p->first * std::min (c, 1.0f - m); // norm = 1
      for (const auto& f : p->second) w[f.first] += t * f.second;
    }
  }
}

// distributed passive aggressive
void train (const std::vector <ex_t> &ex, std::vector <float>& w, const long iter, const float c, const int nf, const int n) {
  w.resize (nf + 1, float (0));
  std::vector <std::thread> thr (n);
  std::vector <std::vector <float> > wv (n, w);
  for (int i = 0; i < iter; ++i) {
    if (n) { // distributed training
      for (int j = 0; j < n; ++j)
        thr[j] = std::thread ([&,i,j] () {
            const size_t nex = ex.size () / n + (ex.size () % n ? 1 : 0);
            if (i) std::copy (std::begin (w), std::end (w), std::begin (wv[j]));
            train1 (&ex[j * nex], &ex[std::min (ex.size (), (j + 1) * nex)], wv[j], c);
          });
      for (int j = 0; j < n; ++j) thr[j].join (); // parallel do
      std::fill (std::begin (w), std::end (w), float (0));
      for (size_t k = 0; k < w.size (); ++k) {
        for (const auto& w_ : wv) w[k] += w_[k];
        w[k] /= n;
      }
    } else // non-distributed training
      train1 (&ex[0], &ex[ex.size ()], w, c);
    std::fprintf (stderr, ".");
  }
}

void test (const char* test, const std::vector <float>& w, const int nf, result_t* r) {
  FILE* reader = std::fopen (test, "r");
  char line[MAX_LINE_LEN];
  r->pp = r->pn = r->np = r->nn = 0;
  for (x_t x; std::fgets (line, MAX_LINE_LEN, reader) != NULL; x.clear ()) {
    char* p = line;
    const int y = static_cast <int> (std::strtol (p, &p, 10));
    float m = 0;
    while (*++p) {
      const int i = strtoi (p);
      if (i > nf) break;
      m += w[i] * strtof (++p);
    }
    if (m >= 0) if (y > 0) ++r->pp; else ++r->np;
    else        if (y > 0) ++r->pn; else ++r->nn;
  }
  std::fclose (reader);
}

int main (int argc, char** argv) {
  if (argc < 5) ::errx (1, "Usage %s train test c iter n", argv[0]);
  const float c    = std::strtof (argv[3], NULL);
  const long  iter = std::strtol (argv[4], NULL, 10);
  const long  n    = std::strtol (argv[5], NULL, 10);
  //
  std::vector <ex_t>  ex;
  std::vector <float>  w;
  int nf = 0;
  result_t r;
  //
  for (int i = 0; i < 40; ++i) power10[i] = ::powf (0.1f, i);
  // read examples
  std::fprintf (stderr, "read: ");
  auto start = std::chrono::steady_clock::now ();
  read_data (argv[1], ex, &nf);
  auto d = std::chrono::steady_clock::now () - start;
  std::fprintf (stderr, "%.3fs\n", std::chrono::duration <double> (d).count ());
  // distributed training
  std::fprintf (stderr, "train: ");
  start = std::chrono::steady_clock::now ();
  train (ex, w, iter, c, nf, n);
  d = std::chrono::steady_clock::now () - start;
  std::fprintf (stderr, "%.3fs\n", std::chrono::duration <double> (d).count ());
  // testing
  std::fprintf (stderr, "test: ");
  start = std::chrono::steady_clock::now ();
  test (argv[2], w, nf, &r);
  d = std::chrono::steady_clock::now () - start;
  std::fprintf (stderr, "%.3fs\n", std::chrono::duration <double> (d).count ());
  std::fprintf (stderr, "acc. %2.3f%% (pp %d) (pn %d) (np %d) (nn %d)\n",
                (r.pp + r.nn) * 100.0 / (r.pp + r.pn + r.np + r.nn),
                r.pp, r.pn, r.np, r.nn);
  return 0;
}

実行してみる．c は分散数別にテストデータに対して最適化（達成できる最高精度を見る目的）注: ローカルディスクからデータを読む実験に差し替えるかも． ローカルディスクで実行しても実験結果の傾向には変化がなかったため，そのままにしておく．

// kdd2010 (algebra); http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary.html#kdd2010%20%28algebra%29
// # of data: 8,407,752 / 510,302 (testing)
// # of features: 20,216,830 / 20,216,830 (testing)
> wget http://www.csie.ntu.edu.tw/\~cjlin/libsvmtools/datasets/binary/kdda.bz2 -O - 2>/dev/null | bzcat | ruby -pe '$_.sub!(/^0/, "-1")' > kdda
> wget http://www.csie.ntu.edu.tw/\~cjlin/libsvmtools/datasets/binary/kdda.t.bz2 -O - 2>/dev/null | bzcat | ruby -pe '$_.sub!(/^0/, "-1")' > kdda.t

> g++ -v 2>&1 | grep version
gcc version 4.8.0 20130118 (experimental) (GCC) 
> g++ -O2 -std=c++11 ppa1.cc -o ppa1 -pthread       

// Intel(R) Xeon(R) CPU E7- 4830  @ 2.13GHz; 32CPU 512GM RAM
// 1 CPU (no threading)
> time ./ppa1 kdda kdda.t 0.1 10 0 
read: 44.586s
train: ..........12.066s
test: 3.817s
acc. 89.497% (pp 437099) (pn 5746) (np 47852) (nn 19605)
./ppa1 kdda kdda.t 0.1 10 0  58.84s user 2.71s system 99% cpu 1:01.64 total

// 1 CPU
> time ./ppa1 kdda kdda.t 0.1 10 1
read: 43.902s
train: ..........18.610s
test: 3.286s
acc. 89.497% (pp 437099) (pn 5746) (np 47852) (nn 19605)
./ppa1 kdda kdda.t 0.1 10 1  64.03s user 2.89s system 99% cpu 1:07.02 total

// 2 CPU
> time ./ppa1 kdda kdda.t 0.2 10 2
read: 44.306s
train: ..........10.447s
test: 2.935s
acc. 89.480% (pp 437969) (pn 4876) (np 48809) (nn 18648)
./ppa1 kdda kdda.t 0.2 10 2  63.44s user 2.86s system 112% cpu 58.907 total

// 4 CPU
> time ./ppa1 kdda kdda.t 0.5 10 4
read: 43.885s
train: ..........7.230s
test: 2.322s
acc. 89.417% (pp 436614) (pn 6231) (np 47775) (nn 19682)
./ppa1 kdda kdda.t 0.5 10 4  66.49s user 2.74s system 126% cpu 54.651 total

// 8 CPU
> time ./ppa1 kdda kdda.t 0.5 10 8
read: 43.602s
train: ..........7.861s
test: 2.321s
acc. 89.324% (pp 437992) (pn 4853) (np 49626) (nn 17831)
./ppa1 kdda kdda.t 0.5 10 8  73.22s user 3.01s system 138% cpu 55.004 total

// 16 CPU
> time ./ppa1 kdda kdda.t 0.5 10 16
read: 42.683s
train: ..........15.302s
test: 2.326s
acc. 89.226% (pp 438419) (pn 4426) (np 50555) (nn 16902)
./ppa1 kdda kdda.t 0.5 10 16  96.91s user 2.99s system 162% cpu 1:01.53 total

// liblinear 1.93
> time train -c 0.01 kdda m2; time predict kdda.t m2 out
.*
optimization finished, #iter = 11
Objective value = -26420.646522
nSV = 5738586
train -c 0.01 kdda m2  200.54s user 3.96s system 96% cpu 3:31.86 total
Accuracy = 89.6197% (457331/510302)
predict kdda.t m2 out  13.37s user 1.52s system 98% cpu 15.090 total

というわけで，予想通りほとんど速くならなかった．実装がヘボいというよりは，単純に分散化される学習部本体に対して訓練例の読み込みやパラメタのコピーのコストが大き過ぎるのが問題．特に，標準形式からデータを読み込む処理 (read) のコストが大き過ぎる．
では素性ベクトルをディスクから読まず動的に生成すればよいかというと，その場合は今度は素性抽出のコストが問題になる．学習部の繰り返し回数が多ければ（学習部のコストが大きければ）意味があるのでは，と思うかもしれないが，それはオンライン学習の進化*1に逆行する考え方で，むしろ学習コストは今後さらに減る方向にあると考えた方が良い．
さらに悪いことには，Iterative Parameter Mixing では分散数を増やせば増やすほど（理論上）収束が遅くなってしまう．8CPU を使った実験ですでにその傾向が認められる．同じ繰り返し回数で分散数を増やしていくと，パラメタ平均化の効果が生じることもあり，ある程度は精度を維持できるのだけど，じきに精度は落ちていく．繰り返し回数を増やせばよいのだけど，そうすると，何のために分散化しているのかよく分からなくなる．また，分散数を増やせば増やすほど，学習部のコストが減り，相対的にパラメタのコピーコストが深刻になる．100CPU 使えば 100倍わくわく，といきたいのだけど，現実はそれにはほど遠いようだ．
素性数が少ないタスクのデータセットであれば，少なくとも学習部に関しては線形にスケールするかもしれない．しかし，そもそも素性数が少ない状況では，訓練例数を増やすことに意味が無い場合が多いような気がするので，役に立ちそうな場面が思いつかない．やはり，論文で分散化している構造学習のように学習部のコストが大きい場面で使うべき手法なのだろう．
これで終わるとあまりに惨めなので，std::strtol と std::strtod を泥臭くいじって入力の読み込みを可能な限り高速化してみた (-DUSE_MY_STRTOI, -DUSE_MY_STRTOF で有効化)．

> g++ -O2 -DUSE_MY_STRTOI -DUSE_MY_STRTOF -std=c++11 ppa1.cc -o ppa1-- -pthread

// Intel(R) Xeon(R) CPU E7- 4830  @ 2.13GHz; 32CPU 512GM RAM
// 1 CPU (no threading)
> time ./ppa1-- kdda kdda.t 0.1 10 0
read: 13.930s
train: ..........12.073s
test: 0.516s
acc. 89.497% (pp 437099) (pn 5746) (np 47852) (nn 19605)
./ppa1-- kdda kdda.t 0.1 10 0  24.98s user 2.64s system 99% cpu 27.690 total

// 1 CPU
> time ./ppa1-- kdda kdda.t 0.1 10 1
read: 13.446s
train: ..........18.420s
test: 0.900s
acc. 89.497% (pp 437099) (pn 5746) (np 47852) (nn 19605)
./ppa1-- kdda kdda.t 0.1 10 1  31.36s user 2.70s system 99% cpu 34.113 total

// 2 CPU
> time ./ppa1-- kdda kdda.t 0.2 10 2
read: 13.449s
train: ..........10.440s
test: 0.654s
acc. 89.480% (pp 437969) (pn 4876) (np 48809) (nn 18648)
./ppa1-- kdda kdda.t 0.2 10 2  30.64s user 2.74s system 129% cpu 25.850 total

// 4 CPU
> time ./ppa1-- kdda kdda.t 0.5 10 4
read: 13.505s
train: ..........7.316s
test: 0.517s
acc. 89.417% (pp 436614) (pn 6231) (np 47775) (nn 19682)
./ppa1-- kdda kdda.t 0.5 10 4  34.78s user 2.68s system 165% cpu 22.567 total

// 8 CPU
> time ./ppa1-- kdda kdda.t 0.5 10 8
read: 13.458s
train: ..........7.866s
test: 0.518s
acc. 89.324% (pp 437992) (pn 4853) (np 49626) (nn 17831)
./ppa1-- kdda kdda.t 0.5 10 8  41.58s user 2.81s system 192% cpu 23.058 total

// 16 CPU
> time ./ppa1-- kdda kdda.t 0.5 10 16
read: 13.542s
train: ..........15.311s
test: 0.519s
acc. 89.226% (pp 438419) (pn 4426) (np 50555) (nn 16902)

訓練例の読み込み部を改善することで訓練時間は半減した*2．分散化による高速化も20%程度認められた．しかし，虚しい・・・．まあ，訓練例を分けるというのがそもそも筋が良くなくて，素性空間を分けるというのが正解なのだろう．

• A Method for Parallel Online Learning by J. Langford

高速化を目的とするのではなく，訓練例がメモリに載らない場合に使うのなら良いのでは？と思う人がいるかもしれない．しかしその場合でも，無駄にたくさんの計算機を使う前に，訓練例をディスク上に圧縮して読み込み易い形で置いてディスクからオンラインで読むことを試す方が良いと感じる．

• オンライン学習でキャッシュを実装してみた - ny23の日記

*1:素性の信頼度を考慮することで収束を高速化するという研究が継続的に行われている．ただし，自分の手持ちの二値素性のデータセットでは（Soft Confidence-Weighted を除いて）知る限り全ての手法でパラメタを調整することで到達できる最高精度が（平均化 Perceptron や PA-I やに比べて）顕著に劣ることもあり，あまり進化しているという印象はないのだけど．

*2:最初のコードの実装が良くなかっただけでは？と思う人がいるかもしれないが，自分の知る限り，世の中に出回っているオンライン学習の実装は，たいてい最初の std::strtol, std::strtod を使ったコードよりも遅い実装になっている．