Autogradで勾配を計算する

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

Autograd

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

$ 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階微分の勾配をそれぞれ求められています。

f:id:ohke:20190726095051p:plain

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

以前の投稿で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)