ニューラルネットワークを実装する [Part 1 線形回帰]

はじめに

最近機械学習の勉強ができてなかったのですが、知人にとあるチュートリアルをおすすめされたので、自分の学習も兼ねて紹介記事を書いていきます。

Notes on machine learning というメモの How to implement a neural networkというシリーズです。
まずこの記事では、Part 1の線形回帰の部分に取り組みます。

peterroelants.github.io

本記事のゴール

本記事では、線形関数とガウシアンノイズで人工的に作り出した観測値をもとに、元の線形関数の予測 (=線形回帰) を行います。
線形回帰には非常に単純なニューラルネットワークを用い、最急降下法でネットワークの重みを最適化していきます。

予測する線形関数の定義

まずは今回の予測対象とする線形関数 $f(x)$ を定義します。元記事に従い
$\displaystyle f(x) = 2x$
とします。さらに、f(x)の値にガウシアンノイズを加えて、観測値 $t$ を定義します。
$\displaystyle t = f(x) + N(0, 0.2)$

# [0, 1]から20点をサンプリング
x = np.random.uniform(0, 1, 20)

# 予測する関数の定義。 f(x) = x * 2
def f(x):
    return x * 2

# ガウシアンノイズを加えて、観測値 t を生成する
# randn関数は正規分布に従う形でランダムな値を出力する関数
noise_variance = 0.2
noise = np.random.randn(x.shape[0]) * noise_variance 
t = f(x) + noise

上記コードで得られた観測値 $t$ と、 $f(x) = 2x$ をプロットすると次のような画像が得られます。
f:id:kanohk:20171015174023p:plain

グラフを出力する部分のコードは省略してあります。
興味のある方は、GitHubにおいてあるNotebookか、元の記事をご覧ください。

モデルとコスト関数の定義

上で定義した関数 $f(x)$ のモデルを次のように定義します。
$\displaystyle {\bf y} = {\bf x} * w$

このモデルは非常に単純なニューラルネットワークで、可視化すると次のようになります。

f:id:kanohk:20171015174957p:plain — 出典: https://peterroelants.github.io/posts/neural_network_implementation_part01/#Linear-regression

このモデルと実際の観測値のずれを表すコスト関数には、次の二乗誤差を用います。(Nは観測値の数で、コード上では20となります。)
$\displaystyle \xi = \sum_{i=1}^N \| t_i - y_i \|^2$

このコスト関数が最小になるような重み $w$ を見つけるのがゴールとなるので、これから行う最適化は次のように記述できます。
$\DeclareMathOperator*{\argmin}{arg\,min} \displaystyle \argmin_w \xi = \argmin_w \sum_{i=1}^N \| t_i - y_i \|^2$

上記関数をPythonで定義します。

# ニューラルネットワークのモデル y = x * w
def nn(x, w):
    return x * w

# コスト関数
def cost(y, t):
    return ((t - y) ** 2).sum()

また、このコスト関数を、実際にノイズを乗せた観測値 $t$ を使ってプロットすると次のようになります。

ws = np.linspace(0, 4, num=100)
cost_ws = np.vectorize(lambda w: cost(nn(x, w), t))(ws)
# (cost_wsをmatplotlibでプロット)

f:id:kanohk:20171015180303p:plain

最急降下法

いよいよコスト関数を最小化することでモデルを最適化していきます。
最急降下法の詳細は他に譲りますが、何度も反復しながら重み $w$ を最適化してゆくアルゴリズムです。
最急降下法での $k$ 回目のイテレーション時の重み $w$ の値を $w(k)$ と書くとき、最急降下法は次のように重みを更新していきます。
$\displaystyle w(k + 1) = w (k) - \Delta w(k)$

勾配 $\Delta w(k)$ は学習率を $\mu$ とするとき次のように定義されます。
$\displaystyle \Delta w = \mu \frac{\partial \xi}{\partial w}$

連鎖率を用いると
$\displaystyle \frac{\partial \xi_i}{\partial w} = \frac{\partial y_i}{\partial w} \frac{\partial \xi_i}{\partial y_i}$
となり、それぞれの微分は
$\displaystyle \frac{\partial \xi_i}{\partial y_i} = \frac{\partial (t_i - y_i)^2}{\partial y_i} = 2 (y_i - ti)$
$\displaystyle \frac{\partial y_i}{\partial w} = \frac{\partial x_i * w}{\partial w} = x_i$
となります。

これを $\Delta w$ に代入して整理すると、
$\displaystyle \Delta w = \mu \cdot 2x_i (y_i - t_i)$
となります。バッチ学習する場合は、単純にこの和を取り
$\displaystyle \Delta w = \mu \cdot 2 \cdot \sum_{i=1}^N x_i (y_i - t_i)$
とします。

上記をコードに落とすとこうなります。イテレーションを4回反復して重みを最適化してみます。

# 勾配を定義
def gradient(w, x, t):
    return 2 * x * (nn(x, w) - t)

# Δwを定義
def delta_w(w_k, x, t, learning_rate):
    return learning_rate * gradient(w_k, x, t).sum()

# 初期の重み
w = 0.1

# 学習率
learning_rate = 0.1

# 最急降下法による更新
num_of_iterations = 4
w_cost = [(w, cost(nn(x, w), t))]
for i in range(num_of_iterations):
    dw = delta_w(w, x, t, learning_rate)
    w = w - dw
    w_cost.append((w, cost(nn(x, w), t)))
    
for i in range(0, len(w_cost)):
    print('w({}): {:.4f} \t cost: {:.4f}'.format(i, w_cost[i][0], w_cost[i][1]))

# w(0): 0.1000 	 cost: 13.6197
# w(1): 1.5277 	 cost: 1.1239
# w(2): 1.8505 	 cost: 0.4853
# w(3): 1.9234 	 cost: 0.4527
# w(4): 1.9399 	 cost: 0.4510

コストが徐々に小さくなり、重みも最初に定義した $f(x) = 2x$ の2に近づいていく様子が観察できます。

今回はコスト関数が非常にシンプルな形なので、2次元のグラフで可視化することができます。
最初の2回の反復部分を可視化してみると、次のようになります。

plt.plot(ws, cost_ws, 'r-')

for i in range(0, len(w_cost) - 2):
    w1, c1 = w_cost[i]
    w2, c2 = w_cost[i + 1]
    plt.plot(w1, c1, 'bo')
    plt.plot([w1, w2], [c1, c2], 'b-')
    plt.text(w1, c1 + 0.5, '$w({})$'.format(i))
    
plt.xlabel('$w$', fontsize=15)
plt.ylabel('$\\xi$', fontsize=15)
plt.title('Gradient descent updates plotted on cost function')
plt.show()

f:id:kanohk:20171015182621p:plain

学習したモデルの可視化

先ほどは最急降下法のイテレーションを4回だけ回しましたが、次は10回まわしてみます。

# 最急降下法のイテレーションを10回まわす
w = 0
num_of_iterations = 10
for i in range(num_of_iterations):
    dw = delta_w(w, x, t, learning_rate)
    w = w - dw

これで得られた重み $w$ を可視化してみると次のようになります。
f:id:kanohk:20171015182926p:plain

20個の観測値をもとに最適化したモデルで、元の関数 $f(x)$ に非常に近いことが確認できます。
今回のモデルではバイアスを定義していないので、このモデルは $x = 0$ の時には $y = 0$ となります。

GitHub

今回使用したJupyter Notebookは下記においてあります。グラフプロット部分等に興味がある方はご覧ください。
github.com

また元の記事にはより詳しい解説もありますので、そちらをご覧になることもお勧めします。
https://peterroelants.github.io/posts/neural_network_implementation_part01

今日も窓辺でプログラム

外資系企業勤めのエンジニアが勉強した内容をまとめておくブログ