Theanoの使い方 (2) 自動微分（2015/5/18）のつづき。今回は、前回までの知識を使ってTheanoで2クラスロジスティック回帰の実装をしてみたい。

Deep Learning Tutorialの最初に紹介されるアルゴリズムもロジスティック回帰なのだが、こちらはMNISTの数字認識を例にしているため多クラスのロジスティック回帰になっているので実装が少し複雑。また、あとで部品として使えるようにオブジェクト指向を使っているためソースコードがわかりにくい。今回はもっとシンプルな実装をしてみよう。

ロジスティック回帰は以下の記事で過去に実装したことがあるが、今回はTheanoの自動微分の機能を使って（確率的ではない）勾配降下法で解いてみたい。

2クラスロジスティック回帰

2クラスのロジスティック回帰は、 $y = 0$ （負例）または $y = 1$ （正例）を分類するタスク。ロジスティック回帰の仮説は

$h_{\theta} (x) = g(\theta^T x)$

$\displaystyle g(z) = \frac{1}{1 + e^{-z}}$

で表される。ロジスティック回帰のコスト関数として以下の交差エントロピー誤差関数を使う。

$\displaystyle J(\theta) = \frac{1}{m} \sum_{i=1}^m \biggl[-y^{(i)} \log(h_{\theta}(x^{(i)})) - (1-y^{(i)}) \log(1-h_{\theta}(x^{(i)}))\biggr]$

前はこの誤差関数の偏微分を自分で解いてから結果の数式を実装したけど、今回はTheanoの自動微分の機能を使う。

ロジスティック回帰については

も参考になる。

コスト関数の自動微分

上の数式をTheanoのシンボルで表現したのが下のコード片。Xは訓練データを表すシンボル（厳密には共有変数）、yは訓練データの正解ラベルを表すシンボル、thetaはパラメータを表すシンボル。Xとthetaの内積の取り方が逆だと思うかもしれないけど、Xは行がサンプル、列が特徴の全データをまとめた行列なのでこの順番でよい。先の数式は小文字のxで1サンプルだけを表している。

g_thetaがコスト関数をパラメータthetaで微分した数式を表すシンボルである。Theanoの自動微分がすごいのはパラメータthetaがスカラーではなくベクトルであってもまったく問題がないこと。あまりにも簡単に書けてしまうのでこれで本当に正しいのか不安になるレベル。

    # 交差エントロピー誤差関数
    h = T.nnet.sigmoid(T.dot(X, theta))
    cost = (1.0 / m) * T.sum(-y * T.log(h) - (1 - y) * T.log(1 - h))

    # コスト関数の微分
    g_theta = T.grad(cost=cost, wrt=theta)

Theanoによる勾配降下法の実装

Theanoを使うと勾配降下法も下のように非常に直観的に書ける。

    # パラメータ更新式
    learning_rate = 0.001
    updates = [(theta, theta - learning_rate * g_theta)]

    # 訓練用の関数を定義
    train_model = theano.function(inputs=[], outputs=cost, updates=updates)

    # 高度な収束判定はせずにiterations回だけ繰り返す
    iterations = 300000
    for iter in range(iterations):
        current_cost = train_model()
        print iter, current_cost

今まで散々使ってきたtheano.function()はupdatesというオプションを指定できる。これを指定しておくと関数の実行後に指定された更新が自動的に行われる。updatesは（共有変数, 更新式のシンボル）のペアが並ぶリスト。今回は、勾配降下法なのでパラメータthetaの更新式をここに入れている。

関数の戻り値にコスト関数を指定しておくのもミソ。こうしておくと関数を呼び出すたびに現在のパラメータでのコスト（誤差）が返される（パラメータのupdates前か後かは要確認）。この誤差が十分小さくなったかを見て収束判定を行うことができる。今回は、高度な収束判定はせずに単純に30万回パラメータ更新を繰り返しているが、後ほどの多クラスロジスティック回帰の実装では、Early stoppingを使った高度な収束判定の実装も紹介されている。

2クラスロジスティック回帰の実装

プログラム全体は下のようになる。学習に使用したデータex2data1.txtはここにある。Xは行方向にサンプル、列方向に特徴量が並ぶ行列で表されているので注意。バイアス項をパラメータthetaに含めるために0次元目に1を加えている。ここら辺は前の記事（2014/4/15）と同じ。

実行すると誤差がだんだん小さくなり、最終的に下のような決定境界が出力される。前にTheanoを使わずに実装した結果（2014/4/15）と同じになることが確認できる。

f:id:aidiary:20150519212003p:plain

(注) データの各次元を平均0、標準偏差1になるように正規化した方がよいかもしれない。sklearn.preprocessing.normalize()が使える。

Theanoによる確率的勾配降下法の実装

さっき実装したのは全訓練データを使って一回だけパラメータを更新する確率的ではない勾配降下法である。ここでは訓練データを1つランダムに選択してパラメータを1回更新する確率的勾配降下法（Stochastic Gradient Descent: SGD）をTheanoで実装してみたい。スクリプトを以下のようにちょっと書き換えるだけで簡単にできる。

ポイントは

SGDではどのような順番でデータを与えるかによって結果が変わるのでデータをシャッフルする
何番目のデータかを表すシンボルindexを定義
コスト関数は全データの和を取らずに1つのデータのみから計算する
train_model()の引数に何番目の訓練データで学習するかを表す引数を追加
各エポックでシャッフルした全訓練データを1つずつ与えてパラメータ更新
全データを使い切ったら次のエポックへ

時間はかかるが学習率はかなり小さめにしている。大きいと訓練データの順番によって誤差が増えたり、発散したりしてしまう。結果は下のような感じ。5000回だとまだ少し誤差が大きいようだ。もう少し回した方がよいかな？でもうまくいくかは運任せかも（笑）

f:id:aidiary:20150520211811p:plain

Deep Learning Tutorialでは、データを複数個まとめたミニバッチという単位で確率的勾配降下法を実行している。これについてはあとで試してみたい。この3つの方法の挙動の違いは、参考文献に挙げたこのサイトで可視化されている。とてもわかりやすいのでおすすめ。

正則化項の導入

最後に正則化項を導入してみよう。正則化を導入したロジスティック回帰のコスト関数は以下の式で表される。パラメータの二乗和が加わっているだけ。

$\displaystyle J(\theta) = \frac{1}{m} \sum_{i=1}^m \biggl[ -y^{(i)} \log(h_{\theta} (x^{(i)})) - (1-y^{(i)}) \log(1-h_{\theta}(x^{(i)})) \biggr] + \frac{\lambda}{2m} \sum_{j=1}^n \theta_j^2$

コスト関数とその微分は下のようなコードになる。微分は自動的に求まるので前と同じコードのままで問題ない。

    # コスト関数を定義
    lam = 1.0
    h = T.nnet.sigmoid(T.dot(X, theta))
    cost = (1.0 / m) * T.sum(-y * T.log(h) - (1 - y) * T.log(1 - h)) + (lam / (2 * m)) * T.sum(theta ** 2)
    # コスト関数の微分
    g_theta = T.grad(cost=cost, wrt=theta)

学習に使用したデータex2data2.txtはここからダウンロードできる。線形分離できないデータなので、前の記事（2014/4/15）でもやったように元の変数の高次の項を導入することでロジスティック回帰でも非線形分離ができるようにしている。

実行結果は

f:id:aidiary:20150524073310p:plain