Autogradで勾配を計算する

最近、大学院で機械学習の講義に通って勉強してます。今回は講義で知ったAutogradについて紹介します。

Autograd

Autogradは勾配を計算するPythonライブラリです。インプットとしてNumPyの行列を渡せる、バックプロパゲーションで計算できるなどの特徴があります。

github.com

$ pip install autograd

Autogradをインポートして使います。

• autograd.numpyはNumPyのラッパになってます
• 勾配を求める関数が grad (elementwise_grad) です (使い方は後述)

import autograd.numpy as np
from autograd import grad, elementwise_grad

使い方はとても簡単で、勾配を求めたい関数をgradに渡すだけです。以下では自然対数関数の勾配を求めてます。

• f'(x=10) = 1/10が得られていることがわかります

import matplotlib.pyplot as plt

def log(x):
    return np.log(x)

grad_log = grad(log)

print("log(10):", log(10.0))
# log(10): 2.302585092994046

print("log'(10):", grad_log(10.0))
# log'(10): 0.1

得られた導関数に対してさらにgradを計算することで、高階微分もできます。

• 入力としてxをarrayで渡す場合、elementwise_gradを使うことで各点の勾配が計算できます

def f(x):
    return 2 * x**2 - 5 * x + 3

x = np.arange(-4.0, 4.0, 0.01)

grad_f1 = elementwise_grad(f)
grad_f2 = elementwise_grad(grad_f1)

plt.plot(x, f(x), label='f')
plt.plot(x, grad_f1(x), label='grad_f1')
plt.plot(x, grad_f2(x), label='grad_f2')
plt.legend()
plt.show()

-4.0〜+4.0の任意の点で、1階微分と2階微分の勾配をそれぞれ求められています。

ロジスティック回帰への適用

以前の投稿でNumPyでロジスティック回帰を実装しましたが、損失関数の勾配の計算をAutogradで置き換えます。

• NumPyでロジスティック回帰を実装する - け日記

ロジスティック回帰の損失関数の勾配も、Autogradを使って求めることができます。

• 損失関数loss_functionをgradで (1階) 微分 -> grad_loss = grad(loss_function)
  • forループを含んでいてもOK!
• 重みwの更新をgrad_lossで計算 -> w = w_t - grad_loss(w_t, X, t) * alpha

from sklearn.datasets import load_iris

# 損失関数
def loss_function(w, X, t):
    loss = 0

    for i in range(t.shape[0]):
        loss -= t[i]*np.dot(w.T, X[:, i]) - np.log(1 + np.exp(np.dot(w.T, X[:, i])))

    return loss

# 最急降下法で学習
def gradient_descent(X, t, alpha=0.001):
    w_t = np.zeros((3, 1))  # wを0で初期化
    loss_t = loss_function(w_t, X, t)
    
    # 損失関数を1階微分
    grad_loss = grad(loss_function)

    w_history, loss_history = [w_t], [loss_t]

    for i in range(1000):
        # 勾配を得てwを更新
        w = w_t - grad_loss(w_t, X, t) * alpha
        loss_t = loss_function(w_t, X, t)
        
        w_history.append(w)
        loss_history.append(loss_t)
        
        w_t = w
        
    return w_history, loss_history

# irisデータセットのロード
iris = load_iris()

# setosaとversicolorの萼片の長さ・幅のみを使う
X = iris.data[iris.target != 2][:, :2].T
X = np.concatenate([X, np.ones((1, X.shape[1]))])  # バイアス項
t = iris.target[iris.target != 2]

# 学習
w_history, loss_history = gradient_descent(X, t)

損失関数の値をプロットすると、収束に向かっていることがわかります。

重みも安定して変化していることがわかります。

cvxpyを使った凸最適化

Pythonで凸最適化を行うための便利なライブラリcvxpyを使う機会がありましたので、使い方を整理しておきます。

凸最適化

凸最適化 (convex optimization) は、制約条件がある中で目的関数の最大化 (または最小化) を行う最適化問題の1つですが、特に以下の特徴を持ちます。

• 目的関数が2次式 (線形を仮定していない)
• 局所的な最大値 (最小値) が大域的な最大値 (最小値) と一致する

変数の集合 とすると、最大化 (または最小化) したい目的関数 、m個の不等式制約 、n個の等式制約 によって一般化されます。

• max. or min.
  • s.t.
  • s.t.

凸最適化は、ラグランジュ関数L ( はラグランジュ乗数) を立式し、

• 

以下のカルーシュ・キューン・タッカー条件 (KKT条件) を使って方程式を立てることで、解析的に導出できます。

• 
• 
• 
• 
• 

cvxpy

cvxpyは上の凸最適化を簡単に計算してくれるPythonライブラリです。

• ドキュメント: https://www.cvxpy.org/
• GitHub: https://github.com/cvxgrp/cvxpy

pipで簡単にインストールできます。

$ pip install cvxpy

cvxpyを使った実装

それではcvxpyを使って凸最適化を行ってみましょう。ここでは以下の問題で最適化を行います。

• min.
• s.t. ,

この問題を解く実装が↓です。

• 変数はVariableで宣言
  • 第1引数shapeで次元数を渡せるので、ベクトルや行列でもOKです
• 目的関数は、最小化ならMinimize、最大化ならMaximizeで、定義
  • Variableとの計算で使える演算子は、こちらを参照
• Problemで目的関数・制約条件を渡してインスタンスを作り、solveメソッドで最適化
  • 制約条件は複数あるためリストにまとめて渡す
• 最適化後の変数は、valueプロパティで取得 ( となってます)
  • 制約条件から、dual_valueプロパティでラグランジュ乗数も取得できます ( ですので2番目の制約条件は使われなかったことになります)

import numpy as np
import cvxpy as cp  # pip install cvxpy

# 変数
x_1 = cp.Variable()
x_2 = cp.Variable()

# 最小化したい関数 (目的関数)
objective = cp.Minimize(
    cp.power(x_1, 2) + 2 * cp.power(x_2, 2)
)

# 制約条件 (不等式制約)
constraints = [
    x_2 >= - x_1 + 1,
    x_2 >= (x_1 - 1) / 2.0
]

# 問題を定義
problem = cp.Problem(objective, constraints)

# 最適化 (戻り値は最適化後に得られた値=最小値)
result = problem.solve()
print(result)  # 0.6666666666666667

# 最小値の時のx_1とx_2の値
print(x_1.value, x_2.value)  # 0.6666666666666667 0.33333333333333337

# 最小値の時のラグランジュ定数
print(constraints[0].dual_value, constraints[1].dual_value)  # 1.3333333333333337 0.0

上では1番目の制約 () 上に最適点がありましたが、制約式上に無い、つまり制約がなかった場合の目的関数の最小値と一致する場合を見ておきます。

目的関数を以下のように書き換えます。この場合、 となり、制約条件は実質不要となります。

# 最小化したい関数 (目的関数)
objective = cp.Minimize(
    cp.power(x_1 - 3, 2) + 2 * cp.power(x_2 - 3, 2)
)

これで最適化すると、予想通り が得られました。さらにラグランジュ乗数が全て0となっており、制約式上に無いことがわかります。

# 最適化
result = problem.solve()
print(result)  # 8.467107550921043e-48

# 最小値の時のx_1とx_2の値
print(x_1.value, x_2.value)  # 3.0 3.0

# 最小値の時のラグランジュ定数
print(constraints[0].dual_value, constraints[1].dual_value)  # 0.0 0.0

もうちょっと見てみます。最小化式を最大化式に無理やり書き換えてみます。この場合、目的関数は無限まで発散します。

# 最小化から最大化する式に変える
objective = cp.Maximize(
    cp.power(x_1, 2) + 2 * cp.power(x_2, 2)
)

実行すると、以下のエラーとなりました。

DCPError: Problem does not follow DCP rules.

まとめ

今回は凸最適化を行うPythonライブラリcvxpyについて、簡単に紹介しました。

Comet.mlでJupyter Notebookの学習を記録・レポートする

Comet.mlを使ってJupyter NotebookなどのPythonの学習を記録・レポートする方法についてまとめます。

[7/18 追記]
本投稿では、稼働環境としてGoogle Colaboratoryを使ってます。環境によってデフォルトで送られるデータに差がありますので、注意が必要です。

Comet.ml

Jupyter Notebookの実験結果を良い感じに・楽に記録したいケースがよくあります。

• ハイパパラメータやシステム設定などの実験条件をトラックしたい
• トレーニングの経過をビジュアライズしたい

Comet.mlはこうしたことを実現するサービスです。

• workspace配下に、複数のprojectを作成し、各プロジェクトでexperimentをトラックする、という構造になってます
• パブリックなプロジェクトであれば無料で利用できます

www.comet.ml

実装

Keras: スパムメッセージをLSTMで分類する - け日記 の実験をComet.mlでレポートできるようにしてみましょう。

ohke.hateblo.jp

comet_mlのインストールとインポート

comet_mlのパッケージをインストールします。このパッケージをノートブック上でインポートすることで

$ pip install comet_ml

インストールしたcomet_mlパッケージをインポートするのですが、注意すべき点としてtensorflowやKerasより前にインポートする必要があります (後でインポートするとエラーになります) 。

from comet_ml import Experiment  # tensorflowやKerasの前にインポートする

import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix
from keras.preprocessing.text import Tokenizer
from keras.preprocessing.sequence import pad_sequences
from keras.models import Sequential
from keras.layers import Input, LSTM, Dense, Embedding

実験の開始

Experimentクラス (APIはこちら) をインスタンス化すると、実験開始となります。
このExperimentクラスがキモとなります。以降はこのインスタンスを介して実験を記録していきます。APIリファレンスを見るとcomet.mlでできることが概ね把握できます。

# Experimentを作る
experiment = Experiment(
    api_key='xxxxxx',  # ワークスペースに紐づくAPIキー
    project_name='lstm_spam'  # プロジェクト名
)

これでComet.mlのワークスペースに、指定したプロジェクト名で新規作成されます。

プロジェクト (ここではlstm-spam) を見ると、experimentsに追加されていきます。

なおAPIキーはこちらで確認できます。

モデリングと学習・推論

あとは以前の投稿と同様にデータロード、前処理、モデル生成を実装していきます。

# ハイパパラメータ
hp = {
    'epochs': 10,  # 学習エポック数
    'batch_size': 32,  # 学習バッチサイズ数
    'optimizer': 'adam',  # 最適化アルゴリズム
    'max_len': 100,  # 1メッセージの最大単語数 (不足分はパディング)
    'embedding_output_dim': 32,  # Embedding層の出力次元数
    'lstm_hidden_units': 16,  # LTSM層の隠れユニット数
    'lstm_activation': 'tanh',  # LSTM層の活性化関数
}

# データロード
dataset_df = pd.read_csv('./SMSSpamCollection', sep='\t', header=None)
dataset_df.rename({0: 'label', 1: 'text'}, axis=1, inplace=True)
dataset_df['category'] = dataset_df.apply(lambda r: 1 if r['label'] == 'spam' else 0, axis=1)

X_train, X_test, Y_train, Y_test = train_test_split(
    dataset_df[['text']], dataset_df[['category']], 
    test_size=0.2, random_state=0
)

# 入力の前処理
tokenizer = Tokenizer()
tokenizer.fit_on_texts(X_train['text'])
x_train = tokenizer.texts_to_sequences(X_train['text'])
x_test = tokenizer.texts_to_sequences(X_test['text'])
x_train = pad_sequences(x_train, maxlen=hp['max_len'])
x_test = pad_sequences(x_test, maxlen=hp['max_len'])

# 出力の前処理
y_train = Y_train['category'].values
y_test = Y_test['category'].values

# モデル生成
vocabulary_size = len(tokenizer.word_index) + 1  # 学習データの語彙数+1
model = Sequential()
model.add(
    Embedding(input_dim=vocabulary_size, output_dim=hp['embedding_output_dim'])
)
model.add(
    LSTM(hp['lstm_hidden_units'], activation=hp['lstm_activation'], return_sequences=False)
)
model.add(Dense(1, activation='sigmoid'))
model.compile(loss='binary_crossentropy', optimizer=hp['optimizer'], metrics=['accuracy'])

学習・推論では、以下のようにそれぞれtrain・testメソッドを呼び出してwith句でくくることで、結果がComet.mlにアップロードされます。

# 学習
with experiment.train():
  history = model.fit(
      x_train, y_train, batch_size=hp['batch_size'], epochs=hp['epochs'],
      validation_data=(x_test, y_test)
  )

# 推論
with experiment.test():
  y_pred = model.predict_classes(x_test)

するとこんな感じでChartsから学習過程などをグラフで見ることができるようになります。

また最終的なaccuracyなどの数字はMetricsで確認できます。

実験条件のロギング

最後に実験条件を送ります。詳しくはAPIドキュメントを見ていただければと思いますが、log_*メソッドで送ることができます。

# 各種条件を履歴に残す
experiment.log_dataset_hash(x_train)  # 学習データ (ハッシュ値)
experiment.log_parameters(hp)  #　ハイパパラメータ

Hyper parametersで送った値が確認できます。特に指定が無い場合でも、最適化アルゴリズムのパラメータなどはデフォルトで残します。

リソース状況などは自動的に記録されており、System Metricsで確認できます。ボトルネックの洗い出しなどはさすがにできませんが、逼迫してないか・有効活用できているかなどのチェックには十分です。

実験の終了

実験完了後はendでクローズします。

# 実験終了
experiment.end()

使う上でのTips

上述した方法で一通りのことはできるようになりましたが、その他実際に使う上で便利なTipsを紹介します。

任意の情報を残したい (log_other)

メトリックスやハイパパラメータ以外に、実験に付随した情報を残したい場合、log_otherメソッドを使うと良いです。
例えば、False negativeとFalse positiveの件数を残したい場合、以下のように記述できます。

• 第1引数がキー、第2引数が値となります

# FN/FPを送る
cm = confusion_matrix(y_test, y_pred)
experiment.log_other('false_negative', cm[1, 0])
experiment.log_other('false_positive', cm[0, 1])

comet.mlでOtherの項目をそれぞれ残っていることがわかります (train_trainable_paramsが送られているのはデフォルトでしょうか) 。

Jupyter Notebookで結果を見たい (display)

Jupyter Notebookでも学習経過を出力したい場合、matplotlibなどで可視化もできますが、displayメソッドでcomet.mlをiframeで表示することもできます。

コードを残したい (set_code)

set_codeメソッドで、実験に使ったコードを残すこともできます。Jupyter NotebookではInに実行したコードをリストで持ってますので、つなげて出力してます。

# コードを出力
code = ''
for cell_in in In:
  code += cell_in + '\n'

experiment.set_code(code)

comet.mlのCodeを開くとビジュアライズされたコードが表示されます。 api_keyの値が消されてて、地味な気遣いを感じます。

まとめ

Jupyter Notebookの学習をComet.mlを使って記録・レポートできるようになりました。