言語処理100本ノック 2015の挑戦記録です。環境はUbuntu 16.04 LTS ＋ Python 3.5.2 :: Anaconda 4.1.1 (64-bit)です。過去のノックの一覧はこちらからどうぞ。
第8章: 機械学習
本章では，Bo Pang氏とLillian Lee氏が公開しているMovie Review Dataのsentence polarity dataset v1.0を用い，文を肯定的（ポジティブ）もしくは否定的（ネガティブ）に分類するタスク（極性分析）に取り組む．
73. 学習
72で抽出した素性を用いて，ロジスティック回帰モデルを学習せよ．
出来上がったコード：
# coding: utf-8
import codecs
import snowballstemmer
import numpy as np
fname_sentiment = 'sentiment.txt'
fname_features = 'features.txt'
fname_theta = 'theta.npy'
fencoding = 'cp1252' # Windows-1252らしい
learn_alpha = 6.0 # 学習率
learn_count = 1000 # 学習の繰り返し数
stemmer = snowballstemmer.stemmer('english')
# ストップワードのリスト http://xpo6.com/list-of-english-stop-words/ のCSV Formatより
stop_words = (
'a,able,about,across,after,all,almost,also,am,among,an,and,any,are,'
'as,at,be,because,been,but,by,can,cannot,could,dear,did,do,does,'
'either,else,ever,every,for,from,get,got,had,has,have,he,her,hers,'
'him,his,how,however,i,if,in,into,is,it,its,just,least,let,like,'
'likely,may,me,might,most,must,my,neither,no,nor,not,of,off,often,'
'on,only,or,other,our,own,rather,said,say,says,she,should,since,so,'
'some,than,that,the,their,them,then,there,these,they,this,tis,to,too,'
'twas,us,wants,was,we,were,what,when,where,which,while,who,whom,why,'
'will,with,would,yet,you,your').lower().split(',')
def is_stopword(str):
'''文字がストップワードかどうかを返す
大小文字は同一視する
戻り値：
ストップワードならTrue、違う場合はFalse
'''
return str.lower() in stop_words
def hypothesis(data_x, theta):
'''仮説関数
data_xに対して、thetaを使ってdata_yを予測
戻り値：
予測値の行列
'''
return 1.0 / (1.0 + np.exp(-data_x.dot(theta)))
def cost(data_x, theta, data_y):
'''目的関数
data_xに対して予測した結果と正解との差を算出
戻り値：
予測と正解との差
'''
m = data_y.size # データ件数
h = hypothesis(data_x, theta) # data_yの予測値の行列
j = 1 / m * np.sum(-data_y * np.log(h) -
(np.ones(m) - data_y) * np.log(np.ones(m) - h))
return j
def gradient(data_x, theta, data_y):
'''最急降下における勾配の算出
戻り値：
thetaに対する勾配の行列
'''
m = data_y.size # データ件数
h = hypothesis(data_x, theta) # data_yの予測値の行列
grad = 1 / m * (h - data_y).dot(data_x)
return grad
def extract_features(data, dict_features):
'''文章から素性を抽出
文章からdict_featuresに含まれる素性を抽出し、
dict_features['(素性)']の位置を1にした行列を返す。
なお、先頭要素は固定で1。素性に対応しない重み用。
戻り値：
先頭要素と、該当素性の位置+1を1にした行列
'''
data_one_x = np.zeros(len(dict_features) + 1, dtype=np.float64)
data_one_x[0] = 1 # 先頭要素は固定で1、素性に対応しない重み用。
for word in data.split(' '):
# 前後の空白文字除去
word = word.strip()
# ストップワード除去
if is_stopword(word):
continue
# ステミング
word = stemmer.stemWord(word)
# 素性のインデックス取得、行列の該当箇所を1に
try:
data_one_x[dict_features[word]] = 1
except:
pass # dict_featuresにない素性は無視
return data_one_x
def load_dict_features():
'''features.txtを読み込み、素性をインデックスに変換するための辞書を作成
インデックスの値は1ベースで、features.txtにおける行番号と一致する。
戻り値：
素性をインデックスに変換する辞書
'''
with codecs.open(fname_features, 'r', fencoding) as file_in:
return {line.strip(): i for i, line in enumerate(file_in, start=1)}
def create_training_set(sentiments, dict_features):
'''正解データsentimentsから学習対象の行列と、極性ラベルの行列を作成
学習対象の行例の大きさは正解データのレビュー数×(素性数+1)。
列の値は、各レビューに対して該当素性がある場合は1、なければ0になる。
列の素性のインデックスはdict_features['(素性)']で決まる。
先頭の列は常に1で、素性に対応しない重みの学習用。
dict_featuresに存在しない素性は無視。
極性ラベルの行列の大きさはレビュー数×1。
肯定的な内容が1、否定的な内容が0。
戻り値：
学習対象の行列,極性ラベルの行列
'''
# 行列を0で初期化
data_x = np.zeros([len(sentiments), len(dict_features) + 1], dtype=np.float64)
data_y = np.zeros(len(sentiments), dtype=np.float64)
for i, line in enumerate(sentiments):
# 素性抽出
data_x[i] = extract_features(line[3:], dict_features)
# 極性ラベル行列のセット
if line[0:2] == '+1':
data_y[i] = 1
return data_x, data_y
def learn(data_x, data_y, alpha, count):
'''ロジスティック回帰の学習
戻り値：
学習済みのtheta
'''
theta = np.zeros(data_x.shape[1])
c = cost(data_x, theta, data_y)
print('\t学習開始\tcost：{}'.format(c))
for i in range(1, count + 1):
grad = gradient(data_x, theta, data_y)
theta -= alpha * grad
# コストとthetaの最大調整量を算出して経過表示（100回に1回）
if i % 100 == 0:
c = cost(data_x, theta, data_y)
e = np.max(np.absolute(alpha * grad))
print('\t学習中(#{})\tcost：{}\tE:{}'.format(i, c, e))
c = cost(data_x, theta, data_y)
e = np.max(np.absolute(alpha * grad))
print('\t学習完了(#{}) \tcost：{}\tE:{}'.format(i, c, e))
return theta
# 素性辞書の読み込み
dict_features = load_dict_features()
# 学習対象の行列と極性ラベルの行列作成
with codecs.open(fname_sentiment, 'r', fencoding) as file_in:
data_x, data_y = create_training_set(list(file_in), dict_features)
# 学習
print('学習率：{}\t学習繰り返し数：{}'.format(learn_alpha, learn_count))
theta = learn(data_x, data_y, alpha=learn_alpha, count=learn_count)
# 結果を保存
np.save(fname_theta, theta)
実行結果：
学習率：6.0 学習繰り返し数：1000
学習開始 cost：0.6931471805599453
学習中(#100) cost：0.4809284917412944 E:0.006248170735186127
学習中(#200) cost：0.43188679850114775 E:0.003599155234198481
学習中(#300) cost：0.4043113376254009 E:0.002616675715766214
学習中(#400) cost：0.38547454091328076 E:0.0020805226234380772
学習中(#500) cost：0.37135664408713015 E:0.0017952496476012821
学習中(#600) cost：0.36017505743644285 E:0.0015873326040173347
学習中(#700) cost：0.35098616931062043 E:0.0014288227472357999
学習中(#800) cost：0.343231725184532 E:0.0013037591670736948
学習中(#900) cost：0.33655507220582787 E:0.0012023865948793643
学習中(#1000) cost：0.33071511988186225 E:0.0011184180264631118
学習完了(#1000) cost：0.33071511988186225 E:0.0011184180264631118
学習した結果の$ \theta $の行列は「theta.npy」に出力します。ファイルはGitHubにアップしています。
NumPyのインストール
準備として、高速な行列演算ができるNumPyというライブラリをインストールする必要があります。幸いAnacondaでインストールされていたようで、何もせずにそのまま使うことができました。オフィシャルサイトはこちらです。
ベクトル化
NumPyによる高速な行列演算が必要になるのは、機械学習の学習作業で必要になる計算量がとても多く、ベクトル化という手法でロジックを行列演算に置き換えながら実装しないと、学習時間がとても長くなってしまうためです。
例えば問題72で出てきた学習に使う次の式は、素性を3,227個も抽出したので、nが3,227になります。
$$ y = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + ... + \theta_{3227}\, x_{3227} $$
この場合、yを1回計算するのに、3,227回の掛け算と3,227回の足し算が必要です。普通に$ \theta_0 $〜$ \theta_{3227} $と$ x_1 $〜$ x_{3227} $をリストにして、for文で3,227回ループしながらそれぞれを掛け算して$ y $に足し込んでいけば良いのですが、学習中はこの式を何回も計算することになるのでかなりの時間がかかってしまいます。
そこで、$ \theta_0 $〜$ \theta_{3227} $と$ x_1 $〜$ x_{3227} $をリストにするのではなく行列にして、行列の内積を取る形で実装します。行列の内積はそれぞれの行列の要素同士を順次掛け算して足し込んでいくことで求められるので、得られる答えは同じです。for文でリスト同士の掛け算と足し込みを繰り返すよりもNumPyで行列の内積を求める方が高速なので、これによって学習を高速化するわけです。行数や列数が1の行列をベクトルと呼ぶので、行列にして演算することを「ベクトル化」と呼んでいます。
なおベクトル化にあたり1つ注意があります。それは$ \theta $が$ \theta_0 $から$ \theta_{3227} $まで計3,228個あるのに対して、$ x $は$ x_1 $から$ x_{3227} $の3,227個で1つ足りないことです。数を合わせないと内積が取れません。そこで、$ x_0 $を用意して式を次のように変形し、$ x_0 $の値を常に1にする形で実装しています。$ x $は素性の数より1つ多くなるので注意してください。
$$ y = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + ... + \theta_{3227}\, x_{3227} (x_0は常に1)$$
行列にした$ X $と$ \theta $は次のようなイメージになります。値は適当です。
これで$ y $は行列の演算による$ X \theta $で求められるようになります。
NumPyによる行列の操作方法は、英語ですがQuickstart tutorialを見るとだいたい分かります。ネットにも解説記事がたくさんありますので、そちらを参照してください。
ベクトル化についての詳細は、junichiroさんの機械学習を1ヵ月で実践レベルにする #6 (Octave 実践編)や、ken5scalさんの[初心者向け]機械学習におけるベクトル化 入門が分かりやすかったです。
ベクトル化でなぜ高速になるのか
ベクトル化してNumPyに行列として計算させても、中でやる計算は同じはずだから高速化しないのでは？と思われるかもしれませんが、NumPyはC/C++並みの速度で演算できるように行列の要素のメモリ配置が工夫されていて、Pythonでループしながら計算するのに比べるとかなり早くなるそうです。
また行列の演算は、内部では要素同士の独立した計算が多いので、逐次で計算するのではなく並列して計算できます。そのため、使うライブラリによってはGPUによる高速化も期待できます。
GPUは当初はグラフィックス専用のチップで、高度なタスクを逐次で高速にこなそうとするCPUコアとは異なり、GPUコアは画面表示のためのピクセルデータの計算タスクを大量に並列実行することに特化していました。
現在のGPUはグラフィックス専用ではなくなり、行列計算のような計算タスクにも使えるようになりました。最近のGPUは数千コアが搭載されているそうなので、GPUコアがこなせる行列演算のようなタスクであれば、コア数が限られているCPUよりもはるかに大量に並列して処理ができます。そのため処理時間の大幅な短縮が期待できるというわけです。
機械学習関連の実装では、ベクトル化は処理を高速化する上で非常に重要な技術と言えそうです。
仮説関数
ここまで例に使ってきた予測のための関数$ h_\theta (x) $を、仮説関数と呼びます。
$$ yの予測値 = h_\theta (x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + ... + \theta_n x_n (x_0は常に1)$$
今回はこの式をもう少し加工して使います。そのまま使うと問題72でも少し話題にしたように、$ y $の予測値が肯定を示す1より大きくなったり、否定を示す0より小さくなったりした時に、いろいろ扱いにくくなるためです。
問題70では、機械学習は結果を正しく予測できたら褒め、不正解なら叱る形で$ \theta $を調整すると書きましたが、実際には後述する目的関数で正解との差を算出して、その差が小さくなるように$ \theta $を調整します。
正しく肯定を予測した場合、$ y $の正解ラベルは1です。でもこの式のままだと、予測値は1だったり10だったり100だったりすることになって、10や100だと正解を正しく予測しているのに正解ラベル1との差が大きくなってしまいます。そこでシグモイド関数という便利な関数を使います。
シグモイド関数
この関数は、値に何を与えても0〜1の間の値に変換してくれるもので、与えた値の大小関係は変換後も維持されるようになっています。この関数の$ z $の部分に仮説関数$ h_\theta (x) $を入れれば、結果が常に0〜1になってくれます。
$$ g(z) = \frac{1}{1 + e^{-z}} $$
グラフに書くとこんな感じです（Mac標準のGrapherに書かせてみましたが、これは便利ですね。ただ、軸ラベルの指定方法が分からないので、そこだけ加工してます）。
$ z $が0近辺で変動した場合はシグモイド関数の結果が0.5付近で大きく変動するのに対し、$ z $がある程度以上大きくなったり小さくなったりすると、シグモイド関数の結果は1付近や0付近でほとんど差がつかなくなります。
この特性のおかげで、仮説関数の結果が1だったり10だったり100だったりした場合でも、$ z $に入れれば1に近い値になり、正解ラベルとの差が小さくなります。
今回はこのシグモイド関数を使った仮説関数を使いました。hypothesis()がその実装です。$ X $と$ \theta $の行列を受け取って予測値を計算するようになっています。
目的関数
続いて、どれくらい正しく予測できているのかを判断するための目的関数の説明です。この目的関数の結果によって、より正しく予測できるように$ \theta $を調整していきます。
仮説関数で計算した予測値と正解のラベルを比較すれば、予測が正しいかがわかります。差が小さければ予測精度が高く、差が大きければ予測精度が低いことを示すので、学習対象のレビューすべてに対して計算して加算し、最後にレビュー数で割って平均を出すのが基本の形です。
通常は単純に差を加算するのではなく、差を二乗した値を加算します。二乗すれば差の符号も気にする必要がなくなり、また差に対する値も大きくなるので、$ \theta $を調整しやすくなるからのようです。
ただし今回のように結果が1（肯定）か0（否定）のどちらかにしかならない予測では、二乗した差を使うだけでは差が最大で1にしかならず、またこの目的関数が凸関数ではなくなってしまって、目的関数を最小値に導くよう$ \theta $を調整するのが難しくなるそうです。
そこで、$ -log() $の特性を利用して、予測と結果が一致した場合は0、一致しない場合は$ \infty $になるよう、次のような目的関数を使います。$ y $が正解ラベル（肯定なら1、否定なら0）、$ h $が仮説関数による予測値です。
$$ -y\,log(h) - (1 - y) log(1 - h) $$
前半で正解が肯定（$ y = 1 $）だった場合の値を計算し、後半で正解が否定（$ y = 0 $）だった場合の値を計算しています。今回はこれを全レビューに対して計算し、それらを加算してレビュー数で割ったものを目的関数として使いました。
目的関数のベクトル化
目的関数の実装も、行列演算を使うことで高速化しています。
普通に実装すると目的関数を計算するためには、レビュー1つずつに対して予測値$ h $を求め、$ -y\,log(h) - (1 - y) log(1 - h) $の値を計算し、それを全レビューの数だけ繰り返して加算していくことになります。これをループで書くと速度が落ちてしまうためです。
まず、全レビューから特徴抽出した行列$ X $を作ります。行列の縦はレビュー数、横は素性数+1になります。+1なのは、ベクトル化のところに書いたように、仮説関数で行列演算を使うために$ x_0 $に相当する常に1の要素が必要なためです。
今回、レビュー数が10,662件、素性数が3,227ですので、$ X $は次のような感じになります。あくまでもイメージで、値は適当です。
$ \theta $はベクトル化のところにも書いた形です。
これで全レビューに対する予測値$ H $が$ X \theta $一発で求められます。
同じように正解ラベルも行列にします。
これで、$ -y\,log(h) - (1 - y) log(1 - h) $の計算も行列の演算でできるようになり、全レビュー分の計算が1回でできるようになります。この形で実装したのがcost()です。
なお、行列$ X $と$ Y $の生成はcreate_training_set()でやっています。
最急降下法
続いて$ \theta $を調整する部分です。
仮説関数で予測し、目的関数で正解との差を計算し、その値が小さくなるように$ \theta $の調整を繰り返します。今回はこの調整に最急降下法を使いました。
最急降下法は、目的関数における現在の個々の$ \theta $の勾配を調べて、目的関数の値が下がる方向へ個々の$ \theta $を調整する方法です。
以下、$ \theta $ が1つしかない場合のイメージです。
まず現在の$ \theta $における①の位置の勾配を調べます。調べた勾配を図では直線で示しています（この直線の長さは後述する学習率に相当します）。そして、目的関数の値が勾配の直線の通りに下がると仮定して、直性の先の位置になるように$ \theta $を変化させます。これで1回目の調整が終わり、$ \theta $が新しい値になって現在位置が②に移動します。
2回目は②の位置の勾配を調べます。以降、この処理を繰り返して③、④と進み、勾配が平らになる位置へ調整していくことで目的関数を最小化します。
前述の目的関数の説明で、凸関数ではないと困ると書きましたが、凸関数ではない場合、このグラフの曲線の途中に最小地点ではないくぼみができます。そうすると最小地点ではないところにハマってしまって最小値まで$ \theta $を導けなくなってしまうわけです。
勾配は偏微分で求められますが求め方は割愛します（というか私もまだ偏微分を思い出せていなくて自力で求められない^^;）。
最終的に次の式で$ \theta $が調整できます。
$$ 新しい\theta_j = 今の\theta_j - \alpha \frac{1}{m} \sum_{i=1}^m (h_\theta (x^{(i)}) - y^{(i)}) x_j ^{(i)} $$
$ m $はレビューの数、$ h_\theta (x^{(i)}) $は$ i $番目のレビューから特徴抽出した$ x^{(i)} $を使って仮説関数で求めた$ i $番目のレビューの予測値、$ y ^ {(i)} $はi番目のレビューの正解ラベル、$ x_j ^{(i)} $は$ i $番目のレビューから特徴抽出した素性$ j $の値、$ \alpha $は学習率と呼ばれるパラメータです。学習率については後述します。
この最急降下法による新しい$ \theta $の計算も、3,228個の$ \theta $に対してループを繰り返して求めるのではなく、行列の演算で高速化しています。勾配の算出はgradient()で、$ \theta $の調整の繰り返しはlearn()で実装しています。
学習率と繰り返し回数
$ \theta $を調整する際の学習率$ \alpha $は、前述の最急降下法の図における勾配の直線の長さに相当します。
大きくすれば$ \theta $の調整量が大きくなるため目的関数の値も大きく減ります。ただし大きくしすぎると目的関数の最小地点を通り越してしまい、目的関数が大きくなったり、それを繰り返して収束しなくなったりします。逆に小さくしすぎると、目的関数の最小地点に向かって着実には進むものの、処理時間が非常に長くります。この辺は試行錯誤が必要な部分です。
また、この$ \theta $の調整を何回繰り返せば学習が完了するのか？というのも試行錯誤が必要な部分です。目的関数は最初のうちはどんどん小さくなりますが、その速度はどんどん落ちていきます。一般的には$ \theta $の調整量が$ 10^{-3} $くらいまで小さくなるのが終わりの目安とのことでした。今回の実装ではlearn()の中で調整量を計算してEとして表示しています。
何回か試行錯誤してみましたが、今回は学習率を6.0として1,000回繰り返してみました。1,000回目の調整でのEは0.001くらいで、目的関数の値もほとんど変化しなくなっていたので、だいたいこの辺で終わって良さそうです。ちなみに学習にかかった時間は手元のパソコンでは1分ほどでした。
機械学習の勉強のススメ（再掲）
今回の記事を書いていて思いましたが、ここまで複雑な内容になると、自分の理解もおぼつかない状態で簡単に説明しようとするのは苦しいです^^;
問題70の機械学習の勉強のススメにも書きましたが、先に世の中のすばらしい教材や記事をぜひ活用してください。
74本目のノックは以上です。誤りなどありましたら、ご指摘いただけますと幸いです。
実行結果には、100本ノックで用いるコーパス・データで配布されているデータの一部が含まれます。