y_uti のブログ

統計、機械学習、自然言語処理などに興味を持つエンジニアの技術ブログです

Soft Confidence-Weighted Learning の実装

[2017-03-20] アルゴリズムと実装の対応など、全体的に説明を補強しました。

Soft Confidence-Weighted (SCW)[1] を実装して、Iris データセットを分類する様子を観察します。SCW-I と SCW-II のほかに Confidence-Weighted (CW)[2], Adaptive Reguralization of Weight Vectors (AROW)[3] も実装して、それぞれの動作を比較します。実装にあたっては『オンライン機械学習』と各アルゴリズムの論文を参考にしました。
www.kspub.co.jp

参考にした論文は以下のとおりです*1*2
[1] J. Wang, P. Zhao, and S. C. H. Hoi. Exact soft confidence-weighted learning. ICML 2012.
[2] K. Crammer, M. Dredze, and F. Pereira. Exact convex confidence-weighted learning. NIPS 2008.
[3] K. Crammer, A. Kulesza, and M. Dredze. Adaptive regularization of weight vectors. NIPS 2009.

データの準備

前回と同様に、Iris データセットから学習データを作成します。データの作成手順については前回の記事を参照してください。利用する学習データは以下のとおりです*3。左は「がく片」の長さと幅を特徴量として "setosa" と "versicolor" をプロットしたものです。線形分離可能なデータの例として利用します。右は花弁の長さと幅を特徴量として "versicolor" と "virginica" をプロットしたもので、線形分離不可能なデータの例として利用します。いずれも、各特徴量の平均が 0 となるようにデータを変換しています。
f:id:y_uti:20170221225648p:plain:w320 f:id:y_uti:20170222094958p:plain:w320

CW による決定境界の学習

はじめに、CW のアルゴリズムを実装して線形分離可能なデータに適用してみます。CW のアルゴリズムは論文[2] の Figure 1 にあります。以下に引用します。

f:id:y_uti:20170320145413p:plain

K. Crammer, M. Dredze, and F. Pereira. Exact convex confidence-weighted learning. NIPS 2008.

CW の実装は以下のとおりです。論文のアルゴリズムとの対応をコメントで示してあります。 \Sigma は一般の分散共分散行列として、式 (10) にしたがいました。

from scipy.stats import norm

def train_cw(eta, max_iter=10000):
    # Initialize
    mu = np.zeros(n_features)
    sigma = np.eye(n_features)
    phi = norm.ppf(eta)
    psi = 1 + phi**2 / 2
    xi = 1 + phi**2

    models = []
    for t in range(max_iter):
        # Receive a training example and true label
        xt, yt = get_training_sample()
        # Compute v and m (11)
        v = xt.dot(sigma).dot(xt)
        m = yt * mu.dot(xt)
        # phi * sqrt(v) > m の場合のみパラメータを更新する (論文[1] の Algorithm 1 による)
        if max(0, phi * np.sqrt(v) - m) > 0:
            # Compute alpha (14)
            alpha = max(0, (-m * psi + np.sqrt((m * phi**2 / 2) ** 2 + v * phi**2 * xi)) / (v * xi))
            # Compute sqrt(u) (『オンライン機械学習』の式 (6.17) による)
            sqrt_u = 2 * v / (alpha * v * phi + np.sqrt((alpha * v * phi) ** 2 + 4 * v))
            # Compute beta (22)
            beta = alpha * phi / (sqrt_u + v * alpha * phi)
            # Update mu and sigma (10)
            mu = mu + alpha * yt * sigma.dot(xt)
            sigma = sigma - beta * sigma.dot(xt)[:,np.newaxis].dot(xt.dot(sigma)[np.newaxis,:])
            # 更新されたパラメータと学習データに対する正解数を記録する
            accuracy = np.count_nonzero(y * X.dot(mu) > 0)
            models.append({'t': t + 1, 'mu': mu, 'sigma': sigma, 'accuracy': accuracy})

    return models

論文と異なる点は、(a) 入力パラメータの  a を 1 に固定した、(b)  \mathtt{max}(0, \; \phi \sqrt{v} - m) > 0 の場合のみパラメータを更新するようにした、(c)  \sqrt{u} の計算で桁落ちを回避するように式変形した、という三点です。(b) は SCW の論文[1] で言及されている内容ですが、CW でも同じ議論が成立するので実装に追加しました。(c) は『オンライン機械学習』の 6.3.3 項の説明にしたがって式 (6.17) を利用したものです*4

プログラム中の get_training_sample 関数は、データセットからランダムに一つのサンプルを取得する関数です。これは前回の記事で実装したものと同じです。

実装したプログラムを用いて Iris データセットを分類します。今回は  \eta を 0.95 に設定してみましたが、この値には特に根拠はありません。

models = train_cw(0.95)

前回の記事で実装した plot_decision_boundary 関数を用いて、決定境界を描画します*5。図のように正しく分類できていることがわかります。

plot_decision_boundary(models[-1]['mu'])

f:id:y_uti:20170316222545p:plain

重みベクトルが更新されていく様子を確認します。各要素について  \pm \sigma の幅を示すように関数を実装しました*6 \Sigma の値が変わるのは重みベクトルが更新されるときに限られるので、このように簡単に線形分離できるデータでは  \Sigma が小さくならないまま学習が終わってしまうようです。

def plot_weights_with_sigma(models):
    labels = [
        'w0 (' + re.sub(r' \(.*\)', '', iris.feature_names[selected_features[0]]) + ')',
        'w1 (' + re.sub(r' \(.*\)', '', iris.feature_names[selected_features[1]]) + ')',
        'w2 (bias)',
    ]
    plt.figure(figsize=(16, 4))
    for i in range(0, 3):
        ts = [m['t'] for m in models]
        mu = [m['mu'][i] for m in models]
        sigma = np.sqrt([m['sigma'][i, i] for m in models])
        plt.fill_between(ts, mu - sigma, mu + sigma, alpha=0.2)
        plt.plot(ts, mu, marker='None', linestyle='-', label=labels[i])
    plt.title('各反復での重みベクトルの変化', fontproperties=fp14)
    plt.xlabel('反復回数', fontproperties=fp11)
    plt.ylabel('重みベクトルの各要素の値', fontproperties=fp11)
    plt.legend(loc='upper right')
    plt.show()

plot_weights_with_sigma(models)

f:id:y_uti:20170316222857p:plain

重みベクトルに対応して決定境界が動いていく様子も確認します。重みベクトルが更新された時刻ごとに、平均  \mu による決定境界と  N(\mu, \Sigma) から 30 本のベクトルをサンプリングした各々による決定境界を描画します。結果は画像ファイルに出力します。

def plot_decision_boundary_by_samples(model, n=10):
    plot_dataset()
    plt.title('Iris dataset (t = {:05d})'.format(m['t']))
    xmin, xmax = plt.xlim()
    ymin, ymax = plt.ylim()
    grid_x, grid_y = np.meshgrid(np.linspace(xmin, xmax, 100), np.linspace(ymin, ymax, 100))
    mu = model['mu']
    sigma = model['sigma']
    weights = np.random.multivariate_normal(mu, sigma, n)
    for i in range(0, n):
        w = weights[i]
        grid_z = w[0] * grid_x + w[1] * grid_y + w[2]
        plt.contour(grid_x, grid_y, grid_z, colors='gray', levels=[0], alpha=0.5, linewidths=0.5)
    grid_z = mu[0] * grid_x + mu[1] * grid_y + mu[2]
    plt.contour(grid_x, grid_y, grid_z, colors='black', levels=[0], linewidths=2)

for i in range(0, len(models)):
    m = models[i]
    plot_decision_boundary_by_samples(m, 30)
    plt.savefig('figures/fig-{:05d}.png'.format(m['t']))
    plt.clf()
plt.close()

出力された画像ファイルから convert コマンドでアニメーション gif を生成したものが下図になります。反復が進むにつれて決定境界のばらつきも小さくなっていく様子がわかります。なお、今回は最大反復回数を 10,000 回に設定しましたが、重みベクトルが更新されたのは 47 回で、最後に更新が発生したのは t = 1,626 でした。

$ cd figures
$ ln -s $(ls -1 | tail -n 1) last.png
$ convert -layers optimize -loop 0 -delay 10 fig-*.png -delay 300 last.png animation.gif

f:id:y_uti:20170316223627g:plain

線形分離不可能なデータセットに対する CW の学習

CW の弱点として、線形分離不可能なデータセットを上手く学習できない問題が知られています。線形分離不可能な学習データに CW を適用して、このことを確認します。先ほどと同様に、決定境界のアニメーションを示します*7。線形分離可能なデータの場合と異なり、重みベクトルが大きく動き続けて収束しない様子を確認できます。
f:id:y_uti:20170316225927g:plain

なお、ここでは 10,000 回の反復計算が最後まで動いた結果を示しましたが、今回の私の実装では、反復計算の途中で  \Sigma が正定値性を満たさない行列になってしまう場合もありました。この点は深く調査できていないのですが、数値計算の精度による問題かもしれません。

SCW による決定境界の学習

次に、CW の制約を緩和したアルゴリズムである SCW を実装します。SCW には SCW-I, SCW-II という二つのバリエーションがありますが、どちらも CW のプログラムから  \alpha の計算式を変更するだけです。

SCW-I での  \alpha の計算式は論文[1] の Proposition 1 にあります。計算式を以下に引用します。CW との違いは、パラメータ  C との min を取る計算が追加される点だけです。

f:id:y_uti:20170320163658p:plain

J. Wang, P. Zhao, and S. C. H. Hoi. Exact soft confidence-weighted learning. ICML 2012.

CW と SCW-I の実装上の差分は以下のとおりです。

3c3
< def train_cw(eta, max_iter=10000):
---
> def train_scw1(eta, cost, max_iter=10000):
26c26
<             alpha = max(0, (-m * psi + np.sqrt((m * phi**2 / 2) ** 2 + v * phi**2 * xi)) / (v * xi))
---
>             alpha = min(cost, max(0, (-m * psi + np.sqrt((m * phi**2 / 2) ** 2 + v * phi**2 * xi)) / (v * xi)))

CW の場合と同様に、SCW-I による 10,000 回の反復後の決定境界、重みベクトルの推移、決定境界のアニメーションを順に示します。学習時のパラメータは  \eta = 0.95,  C = 1 としました。重みベクトルの更新回数は 462 回でした。
f:id:y_uti:20170316231812p:plain
f:id:y_uti:20170316231826p:plain
f:id:y_uti:20170316231835g:plain

SCW-II についても同様に実験を行います。SCW-II での  \alpha の計算式は論文[1] の Proposition 2 にあります。こちらも計算式を以下に引用します。

f:id:y_uti:20170320165720p:plain

J. Wang, P. Zhao, and S. C. H. Hoi. Exact soft confidence-weighted learning. ICML 2012.

CW と SCW-II の差分は以下のとおりです。計算式のとおり、 n \gamma を用いて  \alpha を計算します。 \beta の計算式は CW と同じなので差分はありません。

3c3
< def train_cw(eta, max_iter=10000):
---
> def train_scw2(eta, cost, max_iter=10000):
26c26,27
<             alpha = max(0, (-m * psi + np.sqrt((m * phi**2 / 2) ** 2 + v * phi**2 * xi)) / (v * xi))
---
>             n = v + 1 / (2 * cost)
>             gamma = phi * np.sqrt((phi * m * v) ** 2 + 4 * n * v * (n + v * phi**2))
>             alpha = max(0, (gamma - 2 * m * n - phi**2 * m * v) / (2 * n**2 + 2 * n * v * phi**2))

実行結果は以下のとおりです。こちらも  \eta = 0.95,  C = 1 として実行しました。重みベクトルの更新回数は 912 回でした。今回のデータについては、SCW-I に比べると更新回数が多いほか、 \Sigma もあまり小さくなっていかないようです。決定境界のアニメーションを見ても、SCW-I よりも広がりをもった状態になっているように見えます。
f:id:y_uti:20170316232819p:plain
f:id:y_uti:20170316232827p:plain
f:id:y_uti:20170316232838g:plain

AROW による決定境界の学習

最後に、AROW についても同様の実験を行います。AROW のアルゴリズムは論文[3] の Figure 1 にあります。以下に引用します。

f:id:y_uti:20170320161331p:plain

K. Crammer, A. Kulesza, and M. Dredze. Adaptive regularization of weight vectors. NIPS 2009.

AROW の実装は以下のとおりです。これも CW や SCW と似た実装になりますが、パラメータ更新条件としてヒンジ損失関数を利用する点、 \alpha だけではなく  \beta の計算式も CW や SCW と異なる点に違いがあります。

def train_arow(r, max_iter=10000):
    # Initialize
    mu = np.zeros(n_features)
    sigma = np.eye(n_features)

    models = []
    for t in range(max_iter):
        # Receive a training example and true label
        xt, yt = get_training_sample()
        # Compute margin and confidence
        m = mu.dot(xt)
        v = xt.dot(sigma).dot(xt)
        # パラメータの更新条件をチェックする
        if m * yt < 1:
            # Update using eqs. (7) and (9)
            beta = 1 / (v + r)
            alpha = max(0, 1 - m * yt) * beta
            mu = mu + alpha * yt * sigma.dot(xt)
            sigma = sigma - beta * sigma.dot(xt)[:,np.newaxis].dot(xt.dot(sigma)[np.newaxis,:])
            # 更新されたパラメータと学習データに対する正解数を記録する
            accuracy = np.count_nonzero(y * X.dot(mu) > 0)
            models.append({'t': t + 1, 'mu': mu, 'sigma': sigma, 'accuracy': accuracy})

    return models

なお、AROW では損失関数の値によらず  \beta > 0 が常に成立するので、パラメータの更新条件のチェックを実装しないと、 \Sigma はどんどん小さくなっていきます。これは論文にも以下のように記載があります。一方、CW や SCW は制約のもとで KL divergence を最小化するように  \mu,  \Sigma を更新するので、現在の  \mu,  \Sigma が制約を満たしていれば、更新条件のチェックを行わなくてもパラメータは変わりません*8

The squared-hinge loss yields a conservative (or passive) update for  \mu in which the mean parameters change only when the margin is too small, and we follow CW learning by enforcing a correspondingly conservative update for the confidence parameter  \Sigma, updating it only when  \mu changes.

K. Crammer, A. Kulesza, and M. Dredze. Adaptive regularization of weight vectors. NIPS 2009.

さて、AROW での実行結果は以下のとおりです。パラメータは  r = 1 としました。重みベクトルの更新回数は 3,108 回で、SCW-I, SCW-II よりも頻繁に更新されるようです。一方で、計算が進むにつれて  \Sigma は小さくなっていく傾向があり、決定境界のばらつきも小さいようです。
f:id:y_uti:20170316233010p:plain
f:id:y_uti:20170316233019p:plain
f:id:y_uti:20170316233030g:plain

*1:SCW の論文[1] は http://icml.cc/2012/papers/ からダウンロードできるのですが、このウェブページに記載されている著者順と論文の PDF ファイルの著者順が異なっています。本記事では PDF ファイルの著者順にしたがいました。

*2:[2017-03-20 修正] 本記事の公開時はウェブページの著者順で記載していましたが、論文の著者順に変更しました。

*3:前回の記事で利用したものと同じです。

*4:[2017-03-20 修正] 本記事の公開時のコードは教科書の式 (6.17) ではなく論文の式をそのまま実装したものでしたので、コードを修正しました。

*5:前回は決定境界に加えてヒンジ損失が 0 になる点も描画しましたが、今回は決定境界のみを描画しています。

*6: \Sigma の対角要素の平方根 \sigma としています。

*7:ただし、線形分離可能な場合に比べて重みベクトルの更新回数が非常に多くなるため、10 回更新されるたびに描画するように間引いて出力しました。本記事のこれ以降の結果も同様です。

*8:実際、CW や SCW では、 \mathtt{max}(0, \; \phi \sqrt{v} - m) = 0 であれば  \alpha = \beta = 0 になることを簡単に確かめられます。