About connecting the dots.

data science related trivial things

勾配ブースティングについてざっくりと説明する

最近xgboostがだいぶ流行っているわけですけど，これはGradient Boosting（勾配ブースティング）の高速なC++実装です．従来使われてたgbtより10倍高速らしいです．そんなxgboostを使うにあたって，はてどういう理屈で動いているものだろうと思っていろいろ文献を読んだのですが，日本語はおろか，英語文献でもそんなに資料がなかったので，ある程度概要を把握するのに結構時間を食いました．

そんなわけで，今回は自分の理解の確認も兼ねて，勾配ブースティングについてざっくりと説明してみようかと思います．とはいえ生理解な部分も結構あるので，マサカリが飛んできそう感が大いにしています．腑に落ちる的な理解を優先しているため，数式は一切出てきません．

勾配ブースティングとは

複数の弱学習器を組み合わせるアンサンブル学習には，いくつかの手法がありますが，ブースティングは逐次的に弱学習器を構築していく手法です．逐次的というのは，弱学習器を1つずつ順番に構築していくという意味です．新しい弱学習器を構築する際に，それまでに構築されたすべての弱学習器の結果を利用します．そのためすべての弱学習器が独立に学習されるバギングと比べると，計算を並列化できず学習に時間がかかります．

ブースティングでは，各ステップごとに弱学習器を構築して損失関数を最小化します．その際に，各学習データの扱いはずっと平等ではありません．各学習データのうち，前のステップで間違って識別されたものへのウェイトを重くして，次のステップで間違ったものをうまく識別できるようにしていきます．

各ステップ内でやることは，ようするに損失関数の最小化問題です．これだけ切り出せば，通常の最適化問題とそれほど大きくは変わりません．最適化問題でよく使われる最急降下法やニュートン法なんかをまとめて，勾配降下法ということができます．勾配ブースティングでやっていることは，各ステップのパラメタ最適化の際に，勾配降下法を用いているというだけのことです．もちろん数学的にはいろいろあるわけですけれども，大枠としてはそれだけです．勾配を求めて学習していく，という形をとるので，損失関数をパラメタ行列で微分してあげるのを繰り返して，所定回数に達したらおしまいです．

このあたり，数式的にしっかり追いたいというのであれば，はじパタの11.4や，統計的学習の基礎の10.10，Introduction to Boosted TreesやFEGさんのKDD Cupまとめあたりを参照してください．

はじめてのパターン認識

はじめてのパターン認識

統計的学習の基礎 ―データマイニング・推論・予測―

統計的学習の基礎 ―データマイニング・推論・予測―

  • 作者: Trevor Hastie,Robert Tibshirani,Jerome Friedman,杉山将,井手剛,神嶌敏弘,栗田多喜夫,前田英作,井尻善久,岩田具治,金森敬文,兼村厚範,烏山昌幸,河原吉伸,木村昭悟,小西嘉典,酒井智弥,鈴木大慈,竹内一郎,玉木徹,出口大輔,冨岡亮太,波部斉,前田新一,持橋大地,山田誠
  • 出版社/メーカー: 共立出版
  • 発売日: 2014/06/25
  • メディア: 単行本
  • この商品を含むブログ (3件) を見る

勾配ブースティングのパラメタ

勾配ブースティングでよく使われるのは，弱識別器に決定木をもちいたGBDT(Gradient Boosting Decision Tree)です．xgboostなんかでも，こちらにパラメタ一覧がまとまっています．GBDTの場合，過学習の制御がパラメタ決定の中心になってきます．

例えば各ステップの学習が後続ステップの学習に影響を及ぼしてしまうため，個々のステップの影響を下げて学習速度をゆっくりにするshrinkageと呼ばれるパラメタがあります．shrinkageのパラメタ\nuは0<\nu<1の値を取ります（xgboostだと，このパラメタはetaで表され，デフォルトは0.3になっています）．ある程度小さいほうが，過学習が抑制されて精度が上がります．過学習に関しては，深層学習によくまとまっています．というのは，ニューラルネットワークは非常に表現力が高い手法であるため，過学習に陥りやすいという特性があります．そこで，いかにして過学習を防ぎながら多層ニューラルの最適値を得るかについて，様々な研究の蓄積があります．

深層学習 (機械学習プロフェッショナルシリーズ)

深層学習 (機械学習プロフェッショナルシリーズ)

また個々の弱学習器である，決定木自体の表現力を制御するパラメタとして，木の深さ（xgboostだとmax_depth）や葉の重みの下限（同じくmin_child_weight），葉の追加による損失減少の下限（gamma）といったものがあります．深さが浅いほど，また下限が大きいほど，当然単純な木になりやすいので過学習の抑制に働きます．

あとは各ステップで決定木の構築に用いるデータの割合というパラメタもあります．学習データから非復元抽出したサブサンプルを用いることで，確率的勾配降下法（Stochastic Gradient Descent: SGD）に近い効果が得られると思われます．当然これも過学習抑制で精度向上につながります．SGDについては最近出版されたオンライン機械学習でも詳しく述べられていて，理解の助けになるかと思います．

オンライン機械学習 (機械学習プロフェッショナルシリーズ)

オンライン機械学習 (機械学習プロフェッショナルシリーズ)

またxgboostにはcolsample_bytreeというパラメタもあり，これは各ステップの決定木ごとに用いる特徴量をサンプリングすることだと思われます．これはランダムフォレストで行われているのと同じで*1，特徴量同士の交互作用を考慮した形のモデリングができるという利点があります．

それ以外にも，木の本数（nrounds）もあります．当然木の本数が多いほうが結果も安定するし精度も上がりますが，その一方で学習に時間がかかってしまいます．先ほども述べたように，ブースティングはパラメタの推定時に前のステップの結果を用いるため，各ステップの推定を同時に行うことができません．そのためステップ数が増えること（＝木の本数が増えること）は，計算時間の増加を招きます．

xgboostパッケージを試してみる

ということで，ようやくですがxgboostを使っていくつかパラメタのシミュレーションをしてみます．データはStackingのときにも使ったバイナリデータです．各パラメタを何段階かで変えてみて，予測率の変化をみてみましょう．

smrmkt.hatenablog.jp

値を変えてみたパラメタ以外のものについては，基本的にデフォルトの値に固定しています．パラメタのデフォルト一覧は以下の通りです*2．

パラメタ 値
nrounds 100*3
eta 0.3
gamma 0
max.depth*4 6
min.child.weight 1
subsumple 1
colsumple.bytree 1

ということで，以下に各パラメタの値を動かしたときに，どの程度正解率が変化するかをまとめてみました．そんなに極端に動くわけではありませんね．

nrounds

nrounds accuracy
1 82.9%
10 85.1%
100 84.8%
1000 83.1%

eta

eta accuracy
0.5 83.9%
0.3 84.8%
0.1 85.2%
0.01 84.7%

gamma

gamma accuracy
0 84.8%
0.1 84.8%
0.3 85.4%
0.5 85.3%

max.depth

max.depth accuracy
3 85.3%
6 84.8%
10 85.3%
20 84.5%

min.child.weight

min.child.weight accuracy
0.1 84.6%
1 84.8%
3 85.2%
10 85.4%

subsumple, colsumple.bytree

なぜかこれらについては，パラメタを変えても実際のモデルに変化がなく，何故なのかがよくわからず...

コードは例のごとくgistにあげているので，ご参考までに．

gist.github.com

*1:ランダムフォレストでは，サンプルと特徴量の両方をサンプリングした木を複数作って，その結果の多数決で識別クラスを選びます

*2:max_delta_stepに関しては，このデータ自体が正例か負例に偏りのあるimbalanceなデータではないため，特に検討の対象とはしません

*3:このパラメタにはデフォルトの値がないため，ここでは100にしています．

*4:xgboostのパラメタの_は，R上の変数として使うときは.で置き換えることができるそうです．

Optimizelyのstats engineによる逐次A/Bテスト

ABテストといえば，だいぶ前に有意とか検定とかそのあたりで，データ系の界隈がいろいろと盛り上がっていたのが記憶に残っているトピックなわけですが，今年の1月にABテストの大手Optimizelyのエンジンがリニューアルされてました．これがなかなか興味深いんで，今回はざっくりとその内容をご紹介します*1．

とりあえず元ネタは以下の記事とテクニカルペーパーになります．

blog.optimizely.com

http://pages.optimizely.com/rs/optimizely/images/stats_engine_technical_paper.pdf

以下の内容は，基本的にはそこに書かれている内容の要約になります．

従来のABテストの問題点

これまでの，いわゆる古典的な統計学に従ったABテストの場合，以下のような問題があります．

  • 想定される差分やサンプルサイズについて，事前に見積もっておかないといけない
  • あらかじめ決めたサンプルサイズに達する前に何度も結果を覗くことで，間違った結果を得てしまうことがある
  • たくさんのバリエーションを一度にテストすることで，誤検出率が上がってしまう

サンプルサイズの事前見積もり

ABテストというのは，ようするに比率の差の検定なわけですけれども，実験的にこれを実施するためのお作法として，事前に想定効果を見積もって，それを検出可能なサンプルサイズを決定します．その上で得られたデータに対して，実際に差が出たかどうかを検定します．これは有名なスチューデントのt検定が，ビールの麦芽汁に酵母液をどれくらい入れればよいのかを決定するための手段として生み出されたように，生産現場において条件を変えて実験する，みたいなものにはとてもよく当てはまります．

しかしながら，Webの世界のABテストのように，サンプルは時間とともにどんどん入ってくるような環境だと，逆に足かせになります．穀物生産の現場ならいざ知らず，WebのUIテストのようなもので，事前に想定効果を事前に見積もるのは至難の技です．もちろんベースラインのパターンのコンバージョン率については，あらかじめログ集計でもしておけばいいわけですけれども，新しいパターンのデザインで3%あがるか，5%あがるかなんてよくわかりませんよね．

結果を何度もみることの影響

そしてABテストツールがあって，時間とともにサンプルが増えてコンバージョン率の折れ線グラフが更新されるというのに，事前に決めたサンプルサイズが溜まるまでそれを見ない，というのも現実的ではないわけです．コンバージョン率には揺らぎがあるため，実験期間中の短い間に5%有意ラインを超えることもよくあります．このときたまたま結果を眺めてたとしたら，差が出たと思って実験を早期に打ち切っちゃう（＝偽陽性の結果となるわけです）なんてパターンがかなりあります．

以下の図は，Optimizelyの記事上にある画像を持ってきたものですが，ベースラインに対して新しい実験パターンのコンバージョンの方がずっと優勢なんですが，95%ラインを越えたり越えなかったり，という推移が見て取れます．ごく初期に95%を越えた段階で効果あり，と思って実験をストップしちゃうなんていうのは割とありがちな事態ではないでしょうか．

http://2nwl40151nra3yah86118nxtto0.wpengine.netdna-cdn.com/wp-content/uploads/2015/01/AA-test.png

誤検出率の向上

またABテストの優位性検定でみている有意性というのは，本当は差がないのに差があると判断しちゃう確率を一定以下（例えば5%）に抑える，という枠組みなわけです．しかしながら，同時に10パターンテストして，そのうち1つが（本当に有意な効果を持っていて）有意だという結果が得られたとします．しかし90%有意のテストだったとすると，10個のうち1個は，本当は効果がないのに有意であるという結果になってしまいます．本来はテスト全体で本当は効果がないのに有意である結果が得られる確率を10%にしたいのに，複数テストを同時に走らせることで，それよりもはるかに高い確率になってしまうわけです．

ソリューション

上記の問題を解決するために，彼らは以下の2つの枠組みを用いた新しいテストエンジンを作ったそうです．

逐次検定ベースのテスト

これまで，古典的統計学の枠組みにそった検定を行っていたことによる問題を，逐次検定を用いることで回避しました．逐次検定は，サンプルが追加されるごとに尤度比を計算して，その尤度比が想定していた閾値を超えた時点で有意とみなす，という手法になります．

実験群と統制群の比率の差分は，the law of the iterated logarithmに従って減衰するということらしく，それをモデルに取り入れた形でスコアを定義しました．数式書くの面倒なので，テクニカルペーパーのp7にある(2)式を参照してください．これにより，時間経過に伴う差分の変動を考慮した形で，第1種の過誤を一定に保ったまま繰り返し検定を行えるようにしました．

偽陽性率ではなく偽発見率

上で説明したように，複数パターンでの同時テストは，偽陽性率の上昇という大きな問題を抱えています．これを解決するために，彼らは検定全体で偽陽性が発見される確率=偽発見率を定義し，ベイズの枠組みを用いて信頼区間を定める形をとりました．信頼区間の上限は，以下のような式で定義されます．\thetaを実験群と統制群の差分，\hat{p}をP(\hat{\theta}_n|\theta=0)，\pi_0を帰無仮説が真である事前確率，iはABテストが行われた回数，Nはサンプル数です．

P(\theta=0|\hat{\theta}_{n,i})=P(\theta=0|\hat{p}_i)=\frac{\pi_0\hat{p}_i}{i/N}

すごくざっくりいうと，上の式で定義されるような，従来の偽陽性率よりも厳しい偽発見率という基準を用いることによって，第1種の過誤を減らす形の方策をとりました．これらの方略によって，従来数十%程度あった第1種の過誤が，一桁%にまで減少したそうです．

そんなわけで，すごいざっくりとしたOptimizelyのstats engineの紹介でした．このあたりの資料をいろいろみた挙句，そこまで頑張るならバンディッドでいいんじゃない？　っていう感想が出てきたことはここだけの話です．

*1:細かい数式まわり，私自身もきちっと理解し切れているわけではないので，そこら辺解説してくれる人がいたら嬉しかったりします．

stackingを試してみた

つい先日，stackingについての以下の記事が話題になっていました．

このあたり，私自身は試したことがなかったので，実際に試してみましたよというお話．

コード

Rでちゃちゃっと書きました．データをk分割して，分割サンプルごとに訓練データでロジスティック回帰予測モデル構築→予測結果を訓練データに追加→RandomForestで予測モデルを構築，までが訓練フェーズ．テストフェーズでは構築したロジスティック回帰とRandomForestを使ってテストデータの分類を行いました．

使ったのは手元にあった2000サンプル，分類クラス数は2クラスで，正例負例がそれぞれ1000サンプルでした．素性ベクトルはすべてバイナリの15個の変数になります．

library(randomForest)

# load data
data = read.delim("data/sample.tsv", sep="\t")

# stacking sample
k = 10 # cross validation split number
result = c()
for (i in 1:k) {
  print(i)
  data.splitted = cv(data, k)
  # construct predict model
  data.train = cv.training(data.splitted, 1)
  model.glm = glm(y~., data=data.train, family=binomial)
  data.train = stacking(data.train, model.glm) # commentout this line when not using stacking
  model.rf = randomForest(y~., data=data.train)
  # predict with test data
  data.test = cv.test(data.splitted, 1)
  data.test = stacking(data.test, model.glm) # commentout this line when not using stacking
  model.rf.predict = predict(model.rf, newdata=data.test, type="class")
  result = rbind(result, score(model.rf.predict, data.test$y))
}

上で説明した処理はすべて関数化して，以下のようにまとめました．また判別制度の評価をちゃんとするには交差検定も必要なので，こちらも関数にしてまとめておきました*1．

# create data for k-fold cross validation
cv = function(d, k) {
  n = sample(nrow(d), nrow(d))
  d.randomized = data[n,] # randomize data
  n.residual = k-nrow(d)%%k
  d.dummy = as.data.frame(matrix(NA, nrow=n.residual, ncol=ncol(d)))
  names(d.dummy) = names(d)
  d.randomized = rbind(d.randomized, d.dummy) # append dummy for residuals
  d.splitted = split(d.randomized, 1:k)
  for (i in 1:k) {
    d.splitted[[i]] = na.omit(d.splitted[[i]])
  }
  d.splitted
}

# training data
cv.training = function(d, k) {
  d.train = as.data.frame(matrix(0, nrow=0, ncol=ncol(d[[1]])))
  names(d.train) = names(d[[1]])
  for (i in 1:length(d)) {
    if (i != k) {
      d.train = rbind(d.train, d[[i]])
    }
  }
  d.train
}

# test data
cv.test = function(d, k) {
  d[[k]]
}

# stacking with glm
stacking = function(d, m) {
  d = cbind(d, predict(m, newdata=d, type="response"))
  names(d)[length(d)] = "stacking"
  d
}

# check
score = function(p, r) {
  s = c(0, 0, 0, 0)
  for (i in 1:length(p)) {
    pi = 2-as.integer(p[[i]])
    ri = 2-as.integer(r[i])
    s[pi*2+ri+1] = s[pi*2+ri+1]+1
  }
  s
}

結果

ということで，10分割の交差検定をした結果を見ましょう．以下のように予測値と実際の結果のマトリックスで表示してみます．

# show results
m = matrix(apply(result, 2, sum), 2, 2)
dimnames(m) = list(c("res$n", "res$p"), c("pred$n", "pred$p"))
print(m)
print(m/nrow(data))

stackingなし

stacking処理を無効にするには，単にstacking()関数を読んでいる部分2箇所をコメントアウトしてあげるだけです．結果は．Precisionは850/(850+149)=83.2%で，Recallは850/(850+134)=86.3%でした．

> print(m)
       res$p res$n
pred$p   850   149
pred$n   134   867
> print(m/nrow(data))
       res$p  res$n
pred$p 0.425 0.0745
pred$n 0.067 0.4335

stackingあり

それに対してstackingはどうでしょう．Precisionは827/(827+167)=83.2%で，Recallは827/(827+152)=84.5%でした．むしろ悪くなってる...orz

> print(m)
       res$p res$n
pred$p   827   167
pred$n   152   854
> print(m/nrow(data))
        res$p  res$n
pred$p 0.4135 0.0835
pred$n 0.0760 0.4270

というわけで，とりあえず試してみましたよというお話でした．コード全体はgistにあげました．よろしければどうぞ．

*1:このあたりをまとめるにあたっては，shakezoさんの RでデータフレームをK分割する - shakezoの日記 を参考にさせていただきました．

CaffeのImageNetで特徴量抽出器を動かすまで

前回でCaffeがインストールできたので，とりあえず今回はImageNetの特徴量抽出器を使うまで．Yahoo! JAPAN Tech blogの記事を参考にやってみたら，ハマりどころがたくさんあったので，そのあたりを共有しましょうの会です．ハマりどころを抜けるのに参考にしたのはこことかここになります．特に後者．


準備

とりあえずモデルとかいろいろ落としておきます．このあたり，あるもんだと思ってると普通にNot Foundとかいわれてハマるのが悲しいです．get_caffe_reference_imagenet_model.sh自体をまずは落としておかないといけないとか...

cd ~/caffe/examples/imagenet/
wget https://raw.githubusercontent.com/sguada/caffe-public/master/models/get_caffe_reference_imagenet_model.sh
chmod u+x get_caffe_reference_imagenet_model.sh
./get_caffe_reference_imagenet_model.sh
cd ~/caffe/data/ilsvrc12/
./get_ilsvrc_aux.sh
cd ~/caffe/
wget http://www.vision.caltech.edu/Image_Datasets/Caltech101/101_ObjectCategories.tar.gz
tar xzvf 101_ObjectCategories.tar.gz

パスも通しておかないといけない．

echo "export PYTHONPATH=/Users/smrmkt/Workspace/caffe/python:$PYTHONPATH" >> ~/.zshrc
echo "export DYLD_FALLBACK_LIBRARY_PATH=/usr/local/cuda/lib:$HOME/.pyenv/versions/anaconda-2.0.1/lib:/usr/local/lib:/usr/lib" >> ~/.zshrc
source ~/.zshrc

ImageNetのモデル定義ファイルも，落としてこないといけない．

cd examples/imagenet
wget https://raw.githubusercontent.com/aybassiouny/wincaffe-cmake/master/examples/imagenet/imagenet_deploy.prototxt
cp imagenet_deploy.prototxt imagenet_feature.prototxt

その上で，活性化関数を通す前のfc6層の値を取り出すため，imagenet_feature.prototxtの定義ファイルを変更します．

vim imagenet_feature.prototxt
# edit line 174 & 186
174   top: "fc6wi" # fc6->fc6wi
175   blobs_lr: 1
176   blobs_lr: 2
177   weight_decay: 1
178   weight_decay: 0
179   inner_product_param {
180     num_output: 4096
181   }
182 }
183 layers {
184   name: "relu6"
185   type: RELU
186   bottom: "fc6wi" # fc6->fc6wi
# :q

以上で準備は終わりです．

ImageNetで特徴量抽出

元記事のPythonファイルとほぼ同様のファイルを作成します．違いは12-13行目が，net.~~からcaffe.~~になった点です．この修正をしないと動きません．5-7行目のパスは適宜自分の環境に置き換えてください．ここではcaffeのルートディレクトリからの実行を想定しています．

#! /usr/bin/env python
# -*- coding: utf-8 -*-
import sys, os, os.path, numpy, caffe

MEAN_FILE = 'python/caffe/imagenet/ilsvrc_2012_mean.npy'
MODEL_FILE = 'examples/imagenet/imagenet_feature.prototxt'
PRETRAINED = 'examples/imagenet/caffe_reference_imagenet_model'
LAYER = 'fc6wi'
INDEX = 4

net = caffe.Classifier(MODEL_FILE, PRETRAINED)
caffe.set_phase_test()
caffe.set_mode_cpu()
net.set_mean('data', numpy.load(MEAN_FILE))
net.set_raw_scale('data', 255)
net.set_channel_swap('data', (2,1,0))

image = caffe.io.load_image(sys.argv[1])
net.predict([ image ])
feat = net.blobs[LAYER].data[INDEX].flatten().tolist()
print(' '.join(map(str, feat)))

あとは，引数に画像ファイルを指定して，これを実行するだけ．

python feature.py 101_ObjectCategories/airplanes/image_0001.jpg > tmp.txt

実行すると，4096次元の変数に変換することができます．

cat tmp.txt | tr ' ' '\n' | wc -l
4096

HBase徹底入門

ClouderaさまよりHBase徹底入門を献本いただいたの*1で，だいぶ遅くなりましたが感想をまとめておきたいと思います．

HBase徹底入門 Hadoopクラスタによる高速データベースの実現

HBase徹底入門 Hadoopクラスタによる高速データベースの実現

私自身は，Hadoopに関してはもう何年か触っていますが，HBaseについてはまだ0.89あたりのころに，軽い検証をやっただけだったりします．その間に，かつては火山と呼ばれていたHBaseの安定性もだいぶ向上し，もう火山ではなくなっているということで，めでたい限りです．私の知っている範囲でも，本番プロダクトにHBaseを使うという例をカジュアルに見聞きするようになり，すばらしいことだなぁと思っています*2．

そんな私からみて，この本はO'Reillyの馬本と比べても，初心者にわかりやすく丁寧に書かれていると思います．そもそも馬本が2012/7邦訳刊行と，だいぶ情報が古びてしまっている現在では，Cloudera Managerによる最新のクラスタ構築法がまとまっていたりして，新しくHBaseクラスタを導入しようとしたときに，非常に手助けになるように思います．

個人的には，HBaseの肝といわれるスキーマ設計について詳細に述べられている6章と，それを具体的な例を丁寧に示してくれている7章がとても参考になりました．定番のリバースしたID+IDとか，リージョン分散のためにハッシュID+IDみたいな定番パターンが複数のユースケースでまとめられていたり，タイムスタンプに降順で並べたい別の値を入れるみたいなのも解説されていてよいです．しかしこういう入り組んだスキーマ定義が基本になっているというのは，HBase自体の辛いというかアレな所かなぁという感想にはなりますね... とはいえ，そういう部分をきちっと述べてくれているという意味でも，この本の価値はあるのかなぁと思いました．

*1:ありがとうございます！

*2:そしてHBase1.0も無事リリースされたようで，おめでとうございます．

実装して理解するオンライン学習器(2) - Confidence-Weighted

前回からだいぶ間が空きましたが，その間なんもやってなかったので，いい加減まとめてエントリにしておきます．本当はSCWまでやってからにしたかったんですが，あきらめてCWだけで...

元ネタは前回と同じくICMLの以下の論文です．

Jialei Wang, Peilin Zhao, and Steven C. Hoi. Exact soft confidence-weighted learning. In Proc. of ICML 2012, pages 121–128, 2012.

Confidence-Weighted

モデル

オンライン学習器なので，線形モデルでかつデータ追加ごとに逐次学習を進めていくというモデルになります．CWの特徴は，各パラメタについて平均だけでなく分散も同時に求める点にあります．分散が小さければ小さいほど，より精度の高いパラメタ推定ができている，という理屈になります．

(\mu_{t+1}, \sum_{t+1})=\rm{argmax}_{\mu, \sum} \it{D}_{KL} (\mathcal{N}(\mu, \sum), \mathcal{N}(\mu_t, \sum_t))

\it{D}_{KL}はKLダイバージェンスですね*1．こんな感じで，新しいデータが与えられるごとに，KLダイバージェンスを最小にするような:\mu_{t+1}, \sum_{t+1}を求めていく形になります．

詳細な式展開は論文に譲りますが，最終的にはもう少しシンプルな形の閉形式*2であらわすことができます．あと，こちらでも更新式について書かれています．

ということで，最終的にはPassiveAggressiveと同じような形での実装が可能になります．

実装

ということで，実装式は以下の通りです．パラメタとして\etaがあるので，この値を変えることで，モデルの精度が多少変わります．

#!/usr/bin/env python
#-*-coding:utf-8-*-

from math import sqrt
import numpy as np
from scipy.stats import norm

class ConfidenceWeighted():
    def __init__(self, feat_dim, eta=0.90):
        self.t = 0
        self.m = np.ones(feat_dim)
        self.s = np.diag([1.0]*feat_dim)
        self.eta = eta
        self.phi = norm.cdf(self.eta)**(-1)
        self.psi = 1.0+(self.phi**2)/2.0
        self.zeta = 1.0+self.phi**2

    def predict(self, feats):
        return np.dot(self.m, feats)

    def update(self, y, feats):
        # parameter calculation
        v = np.dot(np.dot(feats, self.s), feats)
        m = y*(np.dot(self.m, feats))
        part = sqrt((m**2)*(self.phi**4)/4.0+v*(self.phi**2)*self.zeta)
        alpha = max(0.0, 1.0/(v*self.zeta)*(-m*self.psi+part))
        u = 0.25*((-alpha*v*self.phi+sqrt((alpha**2)*(v**2)*(self.phi**2)+4.0*v))**2)
        beta = (alpha*self.phi)/(sqrt(u)+v*alpha*self.phi)
        # update parameters
        self.t += 1
        self.m += alpha*y*np.dot(self.s, feats)
        self.s -= beta*np.dot(np.matrix(np.dot(self.s, feats).T*feats), self.s)
        return 1 if np.dot(self.m, feats) > 0 else 0

検証

前回と同じく，libsvmのテストデータから，a1aの訓練データとテストデータを持ってきて使いました*3．

まずはオンライン学習をさせて行ったときの精度の変化です．どのモデルでもほとんど変わらず，しかも精度も低いですね... これならPAのときのほうが精度が良いというションボリな感じの結果です．

http://f.st-hatena.com/images/fotolife/S/SAM/20150215/20150215141307_original.png

そしてテストデータに対しての予測は，こちらはそれなりに高くて70%強というところでしょうか．こちらもPAのほうが高いという...普通はCWのほうが精度がいいはずなので，どこか間違えてるのかもしれません．いろいろションボリですが，まぁ仕方なし．

http://f.st-hatena.com/images/fotolife/S/SAM/20150215/20150215141308_original.png

そんな感じのCW編でした．次はSCWにいきたいものです．

*1:一言でいうと，KLダイバージェンスは分布間の距離みたいなものです．この値が小さいほど，似た分布であると考えることができます．

*2:簡単な加減乗除や関数だけで表せる形の式のことを指します．この形にできれば計算しやすいってことです．詳細はこちらの説明をどうぞ．

*3:前処理等のためにヘルパークラスをいくつか作ってgithubにあげてあります．

Mac OS X 10.10にCaffeをインストールするまで

メモ代わりに手順まとめておきます．基本は install_caffe_osx10.10.md と
CaffeをOS X 10.10 にインストールした // ichyo.jpを参考に，細かい修正を幾つか，という感じです．マシンはmac mini late 2012(core i7 2.3GHz quad core)です．

CUDA

CUDAとドライバーをインストール．しかしGPUがIntel HD Graphics 4000なのでCUDAが使えないことに，後から気がつく... 手順的にはpkgとかdmg落としてきて，そのまま入れるだけ．

BLAS

MacにはもともとBLASが入っているので，何もする必要なし*1．

Anaconda

Python周りのものをあらかた入れる．

brew install pyenv
pyenv install anaconda-2.0.1
pyenv rehash
sudo pyenv local anaconda-2.0.1
sudo pyenv global anaconda-2.0.1

OpenCV

OpenCV入れるためには，Homebrewをちゃんとupdateしないといけなかった．

brew update
brew tap homebrew/science
brew install opencv

Boost

formula変更

Homebrewのformulaを変更してから入れる．

Boost

1.55固定にするようにformulaを修正．

cd /usr/local
git checkout a252214 /usr/local/Library/Formula/boost.rb
C++の標準ライブラリをlibstdc++にする

下記コマンドで，該当ライブラリのformulaを開く．

for x in snappy leveldb protobuf gflags glog szip boost boost-python lmdb homebrew/science/opencv; do brew edit $x; done

開いたformulaに以下の修正を加える*2．

def install
+    # ADD THE FOLLOWING:
+    ENV.append "CXXFLAGS", "-stdlib=libstdc++"
+    ENV.append "CFLAGS", "-stdlib=libstdc++"
+    ENV.append "LDFLAGS", "-stdlib=libstdc++ -lstdc++"
+    # The following is necessary because libtool likes to strip LDFLAGS:
+    ENV["CXX"] = "/usr/bin/clang++ -stdlib=libstdc++"
Boost.python

1.55固定にするようにformulaを修正．

brew edit boost-python

でファイルを開いて以下の修正を加える．

   homepage "http://www.boost.org"
-  url "https://downloads.sourceforge.net/project/boost/boost/1.56.0/boost_1_56_0.tar.bz2"
-  sha1 "f94bb008900ed5ba1994a1072140590784b9b5df"
+  url 'https://downloads.sourceforge.net/project/boost/boost/1.55.0/boost_1_55_0.tar.bz2'
+  sha1 'cef9a0cc7084b1d639e06cd3bc34e4251524c840'
+  revision 2
   head "https://github.com/boostorg/boost.git"

Boostのインストール

for x in snappy leveldb gflags glog szip lmdb homebrew/science/opencv; do brew uninstall $x; brew install --build-from-source --fresh -vd $x; done
brew uninstall protobuf; brew install --build-from-source --with-python --fresh -vd protobuf
brew uninstall boost
brew uninstall boost-python
brew install --build-from-source --fresh -vd boost boost-python

ここまでで前準備終わり．

Caffe

ようやくCaffe本体のインストールに突入．落としてきたら設定をいくつか修正してmakeします．

git clone https://github.com/BVLC/caffe.git
cd caffe
cp Makefile.config.example Makefile.config

Makefile.config

GPUモードが使えないので，CPU_ONLYのコメントアウトを外す．

# CPU-only switch (uncomment to build without GPU support).
- # CPU_ONLY := 1
+ CPU_ONLY := 1

Makefile

10.9やBLAS_INCLUDEあたりで検索して，該当箇所を以下のように修正．

OSバージョンの修正
-   ifneq ($(findstring 10.9, $(shell sw_vers -productVersion)),)
+   ifneq ($(findstring 10.10, $(shell sw_vers -productVersion)),)
        CXXFLAGS += -stdlib=libstdc++
        LINKFLAGS += -stdlib=libstdc++
    endif
BLASのパス修正
else ifeq ($(OSX), 1)
    # OS X packages atlas as the vecLib framework
-   BLAS_INCLUDE ?= /System/Library/Frameworks/vecLib.framework/Versions/Current/Headers/
+   BLAS_INCLUDE ?= /Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.9.sdk/System/Library/Frameworks/Accelerate.framework/Versions/Current/Frameworks/vecLib.framework/Headers/
    LIBRARIES += cblas
-   LDFLAGS += -framework vecLib
+   LDFLAGS += -framework Accelerate
endif

パスの追加

ビルドの前にこっちをやっておかないと，make testかmake runtestあたりで，DYLD_FALLBACK_LIBRARY_PATHが指定されてないことによるエラーが出るはず．.zshrcにパスを追加して*3再読み込み．

vim ~/.zshrc

# pyenv
export PYENV_ROOT="${HOME}/.pyenv"
if [ -d "${PYENV_ROOT}" ]; then
    export PATH=${PYENV_ROOT}/bin:$PATH
    eval "$(pyenv init -)"
fi

# caffe
export PYTHONPATH=/path/to/caffe/python:$PYTHONPATH
export DYLD_FALLBACK_LIBRARY_PATH=/usr/local/cuda/lib:$HOME/.pyenv/versions/anaconda-2.0.1/lib:/usr/local/lib:/usr/lib

source ~/.zshrc

コンパイルとテスト

make all
make test
make runtest

正常に終われば，以下のような結果が出るはず．

[----------] Global test environment tear-down
[==========] 457 tests from 98 test cases ran. (11489 ms total)
[  PASSED  ] 457 tests.

pycaffe

Pythonから呼ぶので，こちらもmakeして，protobufを入れておく．

make pycaffe
pip install protobuf

Pythonから呼ぶ

エラーが出る

ipythonからcaffeを呼ぼうとすると，以下のようなエラーが出る．

caffe Fatal Python error: PyThreadState_Get: no current thread

で，ipythonj自体が強制終了してしまう．caffe-userフォーラムでのやり取りを参考に，boost-pythonを以下のように入れ直したら直った．

brew uninstall boost-python
brew install --build-from-source --fresh -vd boost-python

これにてインストール完了．

最後に

はまったところ

  • OpenCV入れるためにbrew tap homebrew/scienceしたところ，普通にエラーが出たので，brew doctorして，エラーを解消してからbrew updateかけるはめになって割とだるかった
  • MakefileでBLASのパス修正するところ，BLAS_INCLUDEだけじゃなくて，LDFLAGSも修正しなきゃいけないことに気づかなくて修正してなかったら，make testのところで "clang: error: linker command failed with exit code 1" というエラーがでて，原因がよくわからず結構つまった．ちゃんとエラーを読むと，vecLibが見つからないよエラーだったので，手順を見直して把握，という流れ
  • 自分のマシンのGPUのことを考えてなくて，CPU_ONLYフラグをアンコメントしなかったことで何度もmake runtestでこけるという失態．それも "CUDA driver version is insufficient for CUDA runtime version." とかいわれるので，何回かCUDAを入れ直すはめに
  • 最後にPythonからcaffeを呼ぶところでFatal errorが出て，これまた詰まる．いろいろググったけど，結局フォーラムのやり取りをもとに入れ直すだけでよかった

感想

結局最初から最後まで，だいたい4時間くらいかかりました．だいぶくたびれたので，imageNetとか触るのはまた今度ということで．こっちもちゃんとドキュメント読まないと，動かせるようになるまでちょっとかかりそう．

*1:ただし10.10は10.9までとパスが変わっているので，make時のmakefile.configを修正する必要あり．

*2:追加した行には+，削除する行には-を，行頭に入れています．

*3:もちろんbashなら.bashrcで．