近接勾配法応用編その2 ～L1ノルム正則化項によるスパースなロジスティック回帰～

今回は前々回の記事で書いた近接勾配法の応用例第２段ということで、分類問題などで使われるロジスティック回帰の目的関数にL1ノルムの正則化項をつけて、スパースな解を導いてみます。スパースな解を導出するメリットとしては、クラス分類に寄与しない無駄な属性の係数が自動的に0となり、クラス間で有意差がある属性を見つけ出せ、かつ、分類モデルの複雑度も下げられるので過学習の防止も期待できるというメリットがあります。なお、ロジスティック回帰は多クラス分類にも利用できますが、今回は２クラス分類の例で説明します。

ロジスティック回帰

ロジスティック回帰について詳しいことはたくさん文献が出ているのでそちらを参照いただくとして、ここでは概略と最終的な目的関数を示します。
まず、学習データ集合 $\left\{ {\bf x}_i, y_i \right\} (i = 1 \cdots N)$ があるとします。ここで、 ${\bf x}$ は $d$ 次元のデータベクトル、 $y$ は $y \in \left\{0, 1\right\}$ のラベルです。
ここで、 ${\bf x}$ が与えられて $y = 1$ となる事後確率は
$\displaystyle p(y = 1| {\bf x}) = \frac{p(y=1, {\bf x})}{p(y=1, {\bf x}) + p(y = 0, {\bf x})} = \frac{1}{1 + \exp(-a) }= \sigma(a) \tag{1}$
と表せます。ここで、 $a$ は
$\displaystyle a= \ln{\frac{p(y=1, {\bf x})}{p(y=0, {\bf x})}} \tag{2}$
という式であり、ロジット関数とも呼ばれます。 ${\bf x}$ がクラス $y = 1$ に所属する確率が大きければこのロジット関数は大きな値を取るので事後確率 $\sigma(a)$ は1に近づきますし、逆にクラス $y=0$ に所属する確率が大きければロジット関数は小さな値を取るので $\sigma(a)$ は0（つまり、クラス $y=0$ となる事後確率 $1-\sigma(a)$ は1）に近づきます。
そしてこの $a$ を ${\bf x}$ の線形関数としてモデル化したのがロジスティック回帰です。具体的には、

$\displaystyle a = {\bf w}^T {\bf x} + w_0 = \tilde{{\bf w}}^T \tilde {{\bf x}} \tag{3}$
ただし、 $\tilde{{\bf w}}= \left[ w_0 , {\bf w} \right]$ , $\tilde{{ \bf x}}= \left[1, {\bf x}\right]$ *1

となります。
さて、これまでの内容を踏まえて、学習データ集合に対する尤度関数を考えます。上記の議論より、 ${\bf x}_i$ と $\tilde{{\bf w}}$ が与えられて、そのクラスが $y_i = 1$ である尤度は $\sigma(a) = \sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i)$ です。一方で、 $y_i = 0$ である尤度は $(1-\sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i))$ です。このように考えると尤度関数は、
$\displaystyle p({\bf y}| \tilde{{\bf w}})= \prod_{i=1}^{N} \sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i)^{y_i} (1 - \sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i))^{(1- y_i)} \tag{4}$
となります。後は式（4）を最大化する $\tilde{{\bf w}}$ を求めればよいです。ただ、この式（4）のままだと計算しにくいので、以下のような式（4）の負の対数を取った誤差関数を最小化することを考えます。
$\displaystyle E(\tilde{{\bf w}}) = - \ln { p({\bf y}| \tilde{{\bf w}}) } = -\sum_{i=1}^{N} \left\{ {y_i} \ln{\sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i)} + (1- y_i)\ln{(1 - \sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i))} \right\} \tag{5}$
これがロジスティック回帰の目的関数です。

正則化項によるスパース化と近接勾配法の適用

ここまでで普通のロジスティック回帰の目的関数が得られました。このロジスティック回帰をスパースな解にするためには、L1ノルム正則化項をつけた式
$\displaystyle F( \tilde{{\bf w}})= -\sum_{i=1}^{N} \left\{ {y_i} \ln{\sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i)} + (1- y_i)\ln{(1 - \sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i))} \right\} + \| \tilde{{\bf w}} \|_1 \tag{6}$
を近接勾配法により最小化すればよいです。近接勾配法を解くためには、元の目的関数 $E(\tilde{{\bf w}})$ の勾配と正則化項 $\| \tilde{{\bf w}} \|_1$ のproximal operatorが必要でした。
正則化項のproximal operatorは以前の記事で導出したとおりです。また、 $E(\tilde{{\bf w}})$ の勾配は、
$\displaystyle \nabla E(\tilde{{\bf w}}) = \sum_{i=1}^{N} (\sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i) - y_i) \tilde{{\bf x}}_i \tag{7}$
となります。後は、これらに基づき近接勾配法を適用するだけです。

実装と実験

これまでの内容pyhonコードで実装すると下記のようになります。なお、近接勾配法のアルゴリズム部分についても、前々回記事で記したコードから停止条件部分を修正したので、改めて近接勾配法の関数も書きます（以前の記事から変更していない関数については下記に記していません）。

import numpy as np
from numpy.random import *
import matplotlib.pyplot as plt

#近接勾配法
def proximal_gradient(grad_f, prox, gamma, objective, init_x, tol = 1e-9):
    x = init_x
    result = objective(x)
    while 1:
        x_new = prox(x - gamma * grad_f(x), gamma)
        result_new = objective(x_new)
        #停止条件を説明変数差分のノルムに変更
        #if result_new < tol or (np.abs(result - result_new)/np.abs(result) < tol) == True :
        if np.linalg.norm(x_new - x, 2) < tol:
            break;
        x = x_new
        result = result_new
    return x_new, result_new

#l1ノルム正則化ロジスティック回帰
#X:学習用データベクトルを列ベクトルとして並べた行列
#y:ラベルを並べたベクトル
def l1norm_LogisticRegression(X, y, lam):
    sigma = lambda a : 1.0/(1.0+np.exp(-a))
    p=lambda Z, w:sigma(np.dot(Z.T, w))  # ※1下記に補足あり
    X=concatenate((ones((1,size(X,1))), X), 0) #バイアス項の追加
    objective = lambda w: - np.sum( y * np.log(p(X, w)) +
                        (1 - y) *np.log(( 1 - p(X, w)))) + lam * np.sum(np.abs(w))
    grad_E = lambda w: np.dot(X, p(X, w) - y) #※2下記補足あり
    (u, l, v) = np.linalg.svd(X)
    gamma = 1.0/max(l.real*l.real)
    prox = lambda v, gamma:prox_norm1(v, gamma, lam) #proximal operator
    w_init = np.zeros(X.shape[0])
    (w_hat, result) = proximal_gradient(grad_E, prox, gamma, objective, w_init, 1e-12)
    return w_hat, result

繰り返し構文を避けるため、行列・ベクトル演算で勾配計算などを行っています。これを少し補足しておくと、まずソースコード※1に当たる部分は、 $p({\bf Z}, {\bf w}) = \left[\sigma({\bf w}^T {\bf z}_1) , \cdots, \sigma({\bf w}^T {\bf z}_N) \right]^T$ というベクトルを返す関数として定義しています。そして、※2の部分では、
$\displaystyle{\bf X}(p({\bf X}, {\bf w}) - {\bf y}) = \left[ {\bf x}_1, \cdots, {\bf x}_N \right] \left[\sigma({\bf w}^T {\bf x}_1) - y_1, \cdots, \sigma({\bf w}^T {\bf x}_N) - y_N \right]^T= \sum_{i=1}^{N} {\bf x}_i (\sigma({\bf w}^T {\bf x}_i) - y_i)$
というように、行列・ベクトルの演算で式（7）の勾配計算を行う関数を定義しています。
このL1ノルム正則化つきロジスティック回帰を疑似データで実験してみました。疑似データは以下のようなコードで生成しました。

def main():
    #学習データ
    #最初の２次元だけ分類に寄与する属性
    X_train1=np.random.randn(2,200)
    X_train2=np.random.randn(2,200) + ones((2,200))
    X_train = concatenate((X_train1,X_train2),1)
    #残り123次元は分類に寄与しないランダムな属性
    X_train = concatenate((X_train, np.random.randn(123, 400)), 0)
    y_train1 = ones((200,1))
    y_train2 = zeros((200,1))
    y_train = concatenate((y_train1,y_train2),0).squeeze()

    #テストデータ
    X_test1=np.random.randn(2,200)
    X_test2=np.random.randn(2,200)+ones((2,200))
    X_test = concatenate((X_test1,X_test2),1)
    X_test = concatenate((X_test, np.random.randn(123, 400)), 0)
    y_test1=ones((200,1))
    y_test2=zeros((200,1))
    y_test=concatenate((y_test1,y_test2),0).squeeze()

疑似データの次元は125次元ですが、その内最初の2次元だけクラス分類において意味があり、それ以外の123次元は全くのランダムでノイズとなる属性です。また、y=0のクラスにのみオフセットをはかせている形なので、バイアス成分もクラス分類に寄与します。そのため、今回のアルゴリズムを適用すると、最初の2次元の係数およびバイアス項のみが非0となり、それ以外の係数は0になることが期待されます。
このように生成した疑似データに対し、以下のように $\lambda=0, 1, \cdots, 39$ と変化させて、その時の推定精度と $w$ の非0要素の数を求めてみました。

    non_zero_cnt_arr = np.zeros(40)
    acc_arr = np.zeros(40)
    for lam in range(0, 40):
        (w, r) = l1norm_LogisticRegression(X_train, y_train.squeeze(), lam)
        y_result = np.dot(w[1:126], X_test) + w[0]
        y_result[y_result>=0] = 1
        y_result[y_result<0] = 0
        non_zero_cnt_arr[lam] = np.sum(np.abs(w)>1e-12)
        acc_arr[lam] = np.sum(y_result == y_test)/float(y_test.shape[0])

    fig, ax1 = plt.subplots()
    ax1.plot(acc_arr, "r.-")
    ax1.set_ylabel("accuracy")
    ax2=plt.twinx()
    ax2.plot(non_zero_cnt_arr, "bx-")
    ax2.set_ylabel("non zero element count")
    plt.legend()
    plt.show()

その結果は下記のようになりました（赤線が推定精度、青線が係数の非0要素数です）。
f:id:YamagenSakam:20180603015719p:plain

$\lambda = 0$ のときは、すべての要素が非0になりノイズである属性も判別に用いてしまっているので精度が60.5%とあまり良くありません。そこから $\lambda$ を大きくしていくと非0要素数は減少していくことがわかります。一方で推定精度については、適切な $\lambda$ であれば精度が上がっていますが、逆に $\lambda$ が大きすぎると制約がきつくなりすぎるためか精度は下がっていきます。面白いのは、非0要素数が3のとき（クラス分類に寄与する２次元 + バイアス項のみが非0のとき）、推定精度がピークとなると予想していましたが、予想に反して $\lambda=15$ で非0要素数8のときに推定精度が最大（75.5%）になる結果が得られました*2。これは、適切にクラス分類に寄与する属性を見つけられるまで $\lambda$ を上げると、ロジスティック回帰の制約条件としては厳しくなりすぎてしまい結果、精度の低下を招くのだと考えられます。そのため、クラス分類問題とどの属性に有意差があるか見つける問題とは一緒くたにして解くのではなく、分けて解くべきなのだと思います。

まとめ

今回はロジスティック回帰にL1ノルム正則化項をつけることによりスパースな解が得られ、分類器としての精度も向上させられることを見てきました。また、分類に寄与する属性もこのL1ノルム正則化を使うことで見つけられますが、だからと言って制約条件が厳しすぎる $\lambda$ を用いると分類精度低下を招くので注意が必要です。一旦、変数選択の目的でL1ノルム正則化項付きのロジスティック回帰を行って、そこで係数が非0となった要素だけを使って今度は制約なしのロジスティック回帰を再学習させる方法がよいのかも知れません。なお、今回はL1ノルムの正則化を試しましたが、もちろんL1ノルム以外にもグループL1ノルムやL1/L2の混合ノルムも適用できます。
今後の予定としては、もうひとつ近接勾配法の応用として行列データの分類について紹介したり、Alternating Direction Method of Multipliers (ADMM)を調べて書いたり、ICML2018の論文（そろそろ公開？）を読んだりなんかをしたいと思います。

*1:バイアス項 $w_0$ を ${\bf w}$ に統合し、 $d + 1$ 次元のベクトルとして扱います。