Python
online-learning

オンライン分類器の比較

More than 1 year has passed since last update.

動機

前回書いた通り、会社内にデータは全く貯められていない状態です。ですが、将来ログをまともに取得した場合のデータは膨大になることが想定されました。そこで、(時間/空間)計算量を考慮するとオンライン学習アルゴリズムを使うのが最良と判断しました。
(以前のpostも想定しての話を書いています。いろんな意味で残念ですね...orz)
今までオンライン分類器をまともに使った事がなかったため、性能評価も兼ねていくつかの分類器を試してみたというわけです(随分前にですが...)。

オンライン分類器の概要

線形分類器は大体
$w^*:=argmin_wΣ_iL(x^{(i)},y^{(i)},w)+CR(w)$
$L(x^{(i)},y^{(i)},w)$:ロス関数, $R(w)$:正規化項
で表すことができると思います。
オンライン学習では、「データを1つ受け取るたびに逐次的にウェイトを更新する」というのがポイントです。
逐次的にウェイトを更新してデータを捨てる、という事を繰り返すため以下の利点があげられます。

  • 省メモリ
  • 再学習が容易

概要についてはオンライン学習(東大の講義資料?)が参考になると思います。

オンライン分類器といっても、ノイズに強く収束の早いアルゴリズム、自動正規化ができるもの、パラメタの自動調整ができるものなど様々あります。
今回はExact Soft Confidence-Weight Learning(SCW),パラメタ自動調整付きStochastic Gradient Descent(SGD),自動正規化付きSGD,Naive Bayes(NB)を比較してみました。

参考文献

  • No More Pesky Learning Rates, Schaul+, ‘13[pdf]
  • Normalized online learning, Ross+, '13[pdf]
  • Exact Soft Confidence-Weighted Learning, Wang, ICML‘12[pdf]
  • メリーランド大学の講義資料(?)[pdf] ← log lossの偏微分の確認するのに良い資料かと思います

参考コード

以下、作成したコードを記載しました。
多少工夫した点としては、(1)簡単なfeature hashingが出来るようにしたこと、(2)SCWで共分散行列の対角成分のみを保存するようにしたことです。

パラメタ自動調整付きSGD

# coding: utf-8

import numpy as np
import math

class SgdTrainWithAutoLR(object):
    def __init__(self, fname, feat_dim, loss_type):
        self.fname = fname # input file name
        self.weight = None # features weight
        self.feat_dim = 2**feat_dim # max size of feature vector 
        self.bitmask = 2**feat_dim - 1 # mapping dimension
        self.loss_type = loss_type # type of loss function
        self.lambd = None # regularization
        self.gamma = None # learning rate
        self.x = np.zeros(self.feat_dim)
        self.t = 1 # update times
        self.gradbar = np.zeros(self.feat_dim)
        self.grad2bar = np.zeros(self.feat_dim)
        self.tau = np.ones(self.feat_dim)*5 
        self.gbar = np.zeros(self.feat_dim) 
        self.vbar = np.zeros(self.feat_dim) 
        self.hbar = np.zeros(self.feat_dim) 
        self.epsilon = 10**(-9)

        self.loss_weight = 10 # 不均衡率調整用のパラメタ 

        self.update_t = 1

    def train(self,gamma,lambd):
        self.gamma = gamma
        self.lambd = lambd
        self.initialize()
        with open(self.fname,'r') as trainf:
            for line in trainf:
                y = line.strip().split(' ')[0]
                self.features = self.get_features(line.strip().split(' ')[1:])
                y = int(-1) if int(y)<=0 else int(1)  # posi=1, nega=-1にといういつする

                # 数値の予測
                pred = self.predict(self.weight,self.features)                
                grad = y*self.calc_dloss(y*pred)*self.features

                self.gradbar = grad
                self.grad2bar = grad**2

                # update weight
                self.update(pred,y)

                self.t += 1
        print self.weight
        print self.t
        return self.weight

    def initialize(self):
        self.weight = np.zeros(self.feat_dim)

    def get_features(self,data):
        features = np.zeros(self.feat_dim)
        for kv in data:
            k, v = kv.strip().split(':')
            features[int(k)&self.bitmask] += float(v)
        return features

    def predict(self,w,features): #margin
        return np.dot(w,features)

    def calc_loss(self,m): # m=py=wxy
        if self.loss_type == 'hinge':
            return max(0,1-m)          
        elif self.loss_type == 'log':
            if m<=-700: m=-700
            return math.log(1+math.exp(-m))

    # gradient of loss function
    def calc_dloss(self,m): # m=py=wxy
        if self.loss_type == 'hinge':
            res = -1.0 if (1-m)>0 else 0.0 # lossが0を超えていなければloss=0.そうでなければ-mの微分で-1になる
            return res
        elif self.loss_type == 'log':
            if m < 0.0:
                return float(-1.0) / (math.exp(m) + 1.0) # yx-e^(-m)/(1+e^(-m))*yx
            else:
                ez = float( math.exp(-m) )
                return -ez / (ez + 1.0) # -yx+1/(1+e^(-m))*yx

    def update(self, pred, y):
        m = y*pred

        self.gbar *= (1 - self.tau**(-1))
        self.gbar += self.tau**(-1)*self.gradbar

        self.vbar *= (1 - self.tau**(-1))
        self.vbar += self.tau**(-1)*self.grad2bar + self.epsilon

        self.hbar *= (1 - self.tau**(-1))
        self.hbar += self.tau**(-1)*2*self.grad2bar + self.epsilon #

        tmp = self.gbar**2/self.vbar

        # update memory size
        self.tau = (1-tmp)*self.tau+1 + self.epsilon

        # update learning rate
        eta = tmp/self.hbar

        # update weight
        self.update_weight(y, m, eta)

    def update_weight(self,y,m,eta):
        loss = self.calc_loss(m)
        if loss>0.0: # paの場合
            delta = self.calc_dloss(m)*self.features
            self.weight -= eta*y*delta
            self.update_t += 1

    def save_model(self,ofname,w):
        with open(ofname,'w') as f:
            # loss関数のtypeの書き込み
            f.write(self.loss_type+'\n')

            # weightの書き込み
            weight = [str(x).encode('utf-8') for x in w]
            f.write(' '.join(weight)+'\n')

            # 特徴量の次元の書き込み
            f.write(str(self.feat_dim).encode('utf-8'))

自動正規化付きSGD

# coding: utf-8

import numpy as np
import math

class SgdTrainWithAutoNormarize(object):
    def __init__(self, fname, feat_dim, loss_type):
        self.fname = fname # input file name
        self.weight = None # features weight
        self.feat_dim = 2**feat_dim # max size of feature vector 
        self.bitmask = 2**feat_dim - 1 # mapping dimension
        self.loss_type = loss_type # type of loss function

        self.eta = 10.0 # learning rate
        self.features = None
        self.G = np.zeros(self.feat_dim, dtype=np.float64)
        self.t = 1 # update times

        self.epsilon = 10**(-9)

        self.maxfeature = np.zeros(self.feat_dim, dtype=np.float64)
        self.N = 0

    def train(self,gamma,lambd):
        self.gamma = gamma
        self.lambd = lambd
        self.initialize()
        with open(self.fname,'r') as trainf:
            for line in trainf:
                y = line.strip().split(' ')[0]
                self.features = self.get_features(line.strip().split(' ')[1:])
                y = int(-1) if int(y)<=0 else int(1)  # posi=1, nega=-1にといういつする

                # 数値の予測
                pred = self.predict(self.weight,self.features)                

                # weightの正規化
                self.NAG(y,y*pred)

                self.t += 1
        print self.weight
        return self.weight

    def initialize(self):
        self.weight = np.zeros(self.feat_dim,dtype=np.float64)

    def get_features(self,data):
        features = np.zeros(self.feat_dim,dtype=np.float64)
        for kv in data:
            k, v = kv.strip().split(':')
            features[int(k)&self.bitmask] += float(v)
        return features

    def NG(self,y,m):
        """ weightの正規化 """
        idx = np.where( (np.abs(self.features)-self.maxfeature)>0 )

        # weightの調整を行う
        self.weight[idx] *= self.maxfeature[idx]**2/(np.abs(self.features[idx])+self.epsilon)**2

        # max値の更新
        self.maxfeature[idx] = np.abs( self.features[idx] )

        # 係数の更新
        self.N += sum(self.features**2/(self.maxfeature+self.epsilon)**2)

        # weightの更新
        loss = self.calc_loss(m)
        if loss>0.0:
            grad = y*self.calc_dloss(m)*self.features
            self.weight -= self.eta*self.t/self.N*1/(self.maxfeature+self.epsilon)**2*grad

    def NAG(self,y,m):
        """ weightの正規化 """
        idx = np.where( (np.abs(self.features)-self.maxfeature)>0 )

        # weightの調整を行う
        self.weight[idx] *= self.maxfeature[idx]/(np.abs(self.features[idx])+self.epsilon)

        # max値の更新
        self.maxfeature[idx] = np.abs( self.features[idx] )

        # 係数の更新
        self.N += sum(self.features**2/(self.maxfeature+self.epsilon)**2)

        # weightの勾配の計算
        grad = y*self.calc_dloss(m)*self.features
        self.G += grad**2
        self.weight -= self.eta*math.sqrt(self.t/self.N)* \
                        1/(self.maxfeature+self.epsilon)/np.sqrt(self.G+self.epsilon)*grad

    def predict(self,w,features):
        ez = np.dot(w,features)
        #return 1/(1+math.exp(-ez)) #logistic関数にはかけない
        return ez

    def calc_loss(self,m): # m=py=wxy
        if self.loss_type == 'hinge':
            return max(0,1-m)          
        elif self.loss_type == 'log':
            if m<=-700: m=-700
            return math.log(1+math.exp(-m))

    # gradient of loss function
    def calc_dloss(self,m): # m=py=wxy
        if self.loss_type == 'hinge':
            res = -1.0 if (1-m)>0 else 0.0 # lossが0を超えていなければloss=0.そうでなければ-mの微分で-1になる
            return res
        elif self.loss_type == 'log':
            if m < 0.0:
                return float(-1.0) / (math.exp(m) + 1.0) # yx-e^(-m)/(1+e^(-m))*yx
            else:
                ez = float( math.exp(-m) )
                return -ez / (ez + 1.0) # -yx+1/(1+e^(-m))*yx

    def update_weight(self,y,m):
        """ weightの更新は正規化後に行う
        """
        loss = self.calc_loss(m)
        if loss>0.0:
            grad = y*self.calc_dloss(m)*self.features
            self.weight -= self.eta*self.t/self.N*1/(self.maxfeature+self.epsilon)**2*grad

    def save_model(self,ofname,w):
        with open(ofname,'w') as f:
            # loss関数のtypeの書き込み
            f.write(self.loss_type+'\n')

            # weightの書き込み
            weight = [str(x).encode('utf-8') for x in w]
            f.write(' '.join(weight)+'\n')

            # 特徴量の次元の書き込み
            f.write(str(self.feat_dim).encode('utf-8'))

SCW

# coding: utf-8

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

class ScwTrain(object):
    def __init__(self, fname, feat_dim, loss_type, eta, C):
        self.fname = fname # input file name
        self.feat_dim = 2**feat_dim # max size of feature vector 
        self.bitmask = 2**feat_dim - 1 # mapping dimension
        self.loss_type = loss_type # type of loss function
        self.lambd = None # regularization
        self.gamma = None # learning rate
        self.t = 1 # update times
        self.features = np.zeros(self.feat_dim,dtype=np.float64)
        self.mean_weight = np.zeros(self.feat_dim,dtype=np.float64)
        self.sigma_weight = np.ones(self.feat_dim,dtype=np.float64)
        self.alpha = 0
        self.beta = 0
        self.sai = None
        self.pusai = None
        self.u = None
        self.v = None
        self.eta = eta
        self.phai = norm.ppf(eta)
        self.C = C

    def train(self):
        with open(self.fname,'r') as trainf:
            ex_num = 0
            count_y = [0.0, 0.0]
            for line in trainf:
                y = line.strip().split(' ')[0]
                features = line.strip().split(' ')[1:]
                y = int(-1) if int(y)<=0 else int(1)

                # yの予測
                pred = self.predict(self.mean_weight,features)

                # vtの計算
                vt = self.calc_v()

                ex_num += 1

                if self.calc_loss(y)>0:
                    # update weight
                    self.update_param(y,y*pred,vt)
                    if y==-1: count_y[0] += 1
                    else: count_y[1] += 1
        print self.mean_weight
        print "data num=", ex_num
        print "update time=", self.t
        print "count_y=", count_y
        return self.mean_weight

    def initialize(self):
        w_init = np.zeros(self.feat_dim, dtype=np.float64)
        return w_init

    def update_param(self,y,m,v):
        nt = v+float(0.5)/self.C
        gmt = self.phai*math.sqrt( (self.phai*m*v)**2+4.0*nt*v*(nt+v*self.phai**2) )

        self.alpha = max( 0.0 , (-(2.0*m*nt+self.phai**2*m*v)+gmt)/(2.0*(nt**2+nt*v*self.phai**2)) )
        u = 0.25*( -self.alpha*v*self.phai+((self.alpha*v*self.phai)**2 + 4.0*v)**0.5 )**2
        self.beta = self.alpha*self.phai/(u**0.5 + v*self.alpha*self.phai)

        self.mean_weight += self.alpha*y*self.sigma_weight*self.features
        self.sigma_weight -= self.beta*(self.sigma_weight*self.features)**2

        self.t += 1

    def predict(self,w,features):
        val = 0.0
        self.features = np.zeros(self.feat_dim,dtype=np.float64)
        if w != None:
            for feature in features:
                k,v = feature.strip().split(':')
                val += w[int(k) & self.bitmask] * float(v)
                self.features[int(k) & self.bitmask] += float(v)
        return val

    def calc_v(self):
        return np.dot(self.sigma_weight, self.features**2)

    def calc_loss(self,y): # m=py=wxy
        """ 損失関数 """
        res = self.phai * math.sqrt( np.dot(self.sigma_weight, self.features**2) ) \
                - y * np.dot(self.features, self.mean_weight) # sigmaは対角成分のみ保存
        return res      

    def save_model(self,ofname,w):
        with open(ofname,'w') as f:
            f.write(self.loss_type+'\n')

            weight = [str(x).encode('utf-8') for x in w]
            f.write(' '.join(weight)+'\n')

            f.write(str(self.feat_dim).encode('utf-8'))

Naive Bayes

ナイーブベイズは単純なので訓練とテスト一発で書いてますw

# coding: utf-8

import numpy as np
import math

class NB(object):
    def __init__(self, fname, feat_dim):
        self.fname = fname # input file name
        self.feat_dim = 2**feat_dim # max size of feature vector 
        self.bitmask = 2**feat_dim - 1 # mapping dimension
        self.t = 1 # update times
        self.t_select = 1 # times of?@select sample
        self.epsilon = 10**(-9)

        self.N=0.0 # 事例の数の総数
        self.Ny = np.zeros(2, dtype=np.float64) # y=-1, y=1のホルダー(総和のカウント)
        self.Nxy = np.zeros((2,2**feat_dim), dtype=np.float64) # yとxベクトルの組み合わせの出現数の総和

        self.propy = None
        self.propxy = None

    def train(self):
        with open(self.fname,'r') as trainf:
            for line in trainf:
                y = line.strip().split(' ')[0]
                self.features = self.get_features(line.strip().split(' ')[1:])
                #y = int(-1) if int(y)<=0 else int(1)
                y = int(-1) if int(y)<=1 else int(1)  # posi=1, nega=-1にといういつする
                # 事例数をカウント
                self.N += 1

                # yの数をカウントする
                self.incliment_y(y)

                # xとyの組み合わせ数のカウント
                self.incliment_xy(y)

            # 予測に使用する確率値の算出
            self.propy = np.log(self.pred_y()+self.epsilon)
            self.propxy = np.log(self.pred_xy()+self.epsilon)

            for i in xrange(len(self.propy)):
                print self.propxy[i]

    def test(self, ifname):
        with open(ifname,'r') as testf:
            ans = np.zeros(4,dtype=np.int64)
            for line in testf:
                y = line.strip().split(' ')[0]
                features = self.get_features(line.strip().split(' ')[1:])
                y = int(-1) if int(y)<=0 else int(1)

                res = self.test_pred(features)

                # 結果の表の計算
                ans = self.get_ans(ans, y, res)
            print ans

    def get_features(self,data):
        features = np.zeros(self.feat_dim)
        for kv in data:
            k, v = kv.strip().split(':')
            features[int(k)&self.bitmask] += float(v)
        return features

    def incliment_y(self,y):
        if y==-1: self.Ny[0] += 1.0
        else: self.Ny[1] += 1.0  

    def incliment_xy(self,y):
        if y==-1:
            self.Nxy[0] += (self.features!=0)*1.0
        else:        
            self.Nxy[1] += (self.features!=0)*1.0

    def test_pred(self,features):
        res = np.zeros(2,dtype=np.float64)
        for i in xrange(len(self.Ny)):
            res[i] = self.propy[i] \
                    + sum( self.propxy[i]*((features!=0)*1.0) )
        if res[0]>res[1]:
            return -1
        else:
            return 1

    def predict(self):
        res = np.zeros(2,dtype=np.float64)
        predy = np.log(self.pred_y()) # yの確率値の算出
        predx = np.log(self.predxy()) # yの条件付きxの確率値の算出

        res = np.zeros(2,dtype=np.float64)
        for i in xrange(len(self.Ny)):
            res[i] = predy[i]+sum(predx[i])
        if res[0]>res[1]:
            return -1
        else:
            return 1

    def pred_y(self):
        return self.Ny/self.N

    def pred_xy(self):
        res = np.zeros((2,self.feat_dim),dtype=np.float64)
        for i in xrange(len(self.Ny)):
            if self.Ny[i]==0.0:
                res[i] = 0.0
            else:
                res[i] = self.Nxy[i]/self.Ny[i]
        return res

    def get_ans(self,ans,y,res):
        if y==1 and res==1: # 真陽性
            ans[0] += 1
        elif y==1 and res==-1: # 偽陰性
            ans[1] += 1
        elif y==-1 and res==1: # 偽陽性
            ans[2] += 1
        else: # 真陰性
            ans[3] += 1

        return ans

if __name__=='__main__':
    trfname = 'training data file name'
    tefname = 'test data file name'

    bf = BF(trfname, 6)
    bf.train()
    bf.test(tefname)

実験結果

実験に使用したデータはLibsvmのdatasetにあるa9aとcovtype.binaryを使いました。
a9aは2値のデータ、covtypeはスケールが6,000倍くらい違うようなデータがボコボコあるデータです。
NBに連続量を入力する場合は単純に1以上の値は1として量子化しました。
また、covtype.binaryはtest setが用意されていないため、40万件のデータをサンプリングし訓練データとして使用しました。
ハイパーパラメタの調整は特にしていません。それから、イテレートはしていません。
結果は以下のとおりです。
result.png

まとめ

変数ごとのスケールが全然違うデータを使う場合は正規化が重要な事がわかります。
あと、SCWは少ないサンプル数で収束するっぽいな、っていうのがわかると思います。
また、オンライン学習はデータの入力順によって結果が変わってきますので、データをシャッフルして実験も行いました。SCWと自動正規化付きのSGDのみ実験しました。
結果としては、SGDは最大3%程度の変動がありましたが比較的accuracyは安定していました。
一方SCWは最大15%の変動があり、データを1回しかなめない場合はデータ順に比較的左右されやすいようです。
そういえば、この辺りはCROSS2014(機械学習CROSS 前半資料)でも話が出ていましたね。

ちなみにSGDのnormalize、lasso、ridgeなどはVowpal Wabbitに全て実装されているのでVowpal Wabbit使用するのをオススメします。計算速度も爆速です。

間違いありましたら、ご指摘いただけますと助かります。