Python3で多変数関数の最急降下法を実装

はじめに

ニューラルネットワークについて勉強中です.
その確認ついでに,Python3で多変数関数の最急降下法を実装していきます.

最急降下法

最急降下法とは

関数の傾きを調べることで,その関数の最小値を探索する手法です.
具体的には,傾きが最も急な方向に向かって降下させていき,変数の値を順次変化させることで,関数の出力を最小値に近づけていきます.

再急降下法のアルゴリズム

具体的には,以下の通りです.

対象とする関数を fとし,関数 fの引数のベクトルを X = (x_{0}, x_{1}, ..., x_{n}) とします.
また, \eta を学習率とします.

まず, fを引数である X = (x_{0}, x_{1}, ..., x_{n}) の各要素で偏微分し,以下のようなベクトルを計算します.
これを,勾配と言います.
 \displaystyle (\frac{\partial f}{\partial x_{0}},  \frac{\partial f}{\partial x_{1}},  ..., 
 \frac{\partial f}{\partial x_{n}})


次に,上で計算した勾配の各要素に学習率を乗じたものを計算します.
 \displaystyle \eta(\frac{\partial f}{\partial x_{0}}),  \eta(\frac{\partial f}{\partial x_{1}}),  ... \eta(\frac{\partial f}{\partial x_{n}})


最後に,全ての  \displaystyle x_{k} (k=0,1,...,n) から,上で計算したものを引き,新たな \displaystyle x_{k}とします.
繰返し回数が \displaystyle t の時の  \displaystyle x_{k} \displaystyle x_{k}^{(t)} と表すとき
  
\displaystyle

\begin{align}
\left\{
\begin{array}{ll}
x_{0}^{(t+1)} =  x_{0}^{(t)} - \eta(\frac{\partial f}{\partial x_{0}^{(t)}}) \\
x_{1}^{(t+1)} =  x_{1}^{(t)} - \eta(\frac{\partial f}{\partial x_{1}^{(t)}}) \\
... \\
x_{n}^{(t+1)} =  x_{n}^{(t)} - \eta(\frac{\partial f}{\partial x_{n}^{(t)}})
\end{array}
\right.
\end{align}
です.

最急降下法の実装

使用するモジュールのインポート

今回は,numpyを使用します.

import numpy as np

勾配を計算する関数の実装

def calc_gradient(f, X):
    """
    calc_gradient
    偏微分を行う関数
    関数fを変数xの各要素で偏微分した結果をベクトルにした勾配を返す
    
    @params
    f: 対象となる関数
    X: 関数fの引数のベクトル(numpy.array)
    
    @return
    gradient: 勾配(numpy.array)
    """
    
    h = 1e-4
    gradient = np.zeros_like(X)
    
    # 各変数についての偏微分を計算する
    for i in range(X.size):
        store_X = X
        
        # f(x+h)
        X[i] += h
        f_x_plus_h = f(X)

        X = store_X
        
        # f(x-h)
        X[i] -= h
        f_x_minus_h = f(X)
        
        X = store_X
        
        # 偏微分
        gradient[i] = (f_x_plus_h - f_x_minus_h) / (2 * h)
        
    return gradient

最急降下法を行う関数

def gradient_descent(f, X, learning_rate, max_iter):
    """
    gradient_descent
    最急降下法を行う関数
    
    @params
    f: 対象となる関数
    X: 関数fの引数のベクトル(numpy.array)
    learning_rate: 学習率
    max_iter: 繰り返し回数
    
    @return
    X: 関数の出力を最小にする(であろう)引数(numpy.array)
    """
    
    for i in range(max_iter):
        X -= (learning_rate * calc_gradient(f, X))
        print("[{:3d}] X = {}, f(X) = {:.7f}".format(i, X, f(X)))
        
    return X

実験

適切な学習率を設定した場合

繰返しを重ねるごとに,f(X)の値が0に近づいていきます

f = lambda X: X[0]**2 + X[1]**2
X = np.array([3.0, 4.0])
gradient_descent(f, X, learning_rate=0.1, max_iter=100)
[  0] X = [ 2.699995  3.599995], f(X) = 20.2499370
[  1] X = [ 2.4299905  3.2399905], f(X) = 16.4023923
[  2] X = [ 2.18698645  2.91598645], f(X) = 13.2858867
[  3] X = [ 1.9682828  2.6243828], f(X) = 10.7615223
[  4] X = [ 1.77144952  2.36193952], f(X) = 8.7167917
[  5] X = [ 1.59429957  2.12574057], f(X) = 7.0605641
[  6] X = [ 1.43486461  1.91316151], f(X) = 5.7190234
[  7] X = [ 1.29137315  1.72184036], f(X) = 4.6323789
...
[ 50] X = [ 0.01386542  0.01850382], f(X) = 0.0005346
...
[ 96] X = [  5.93079900e-05   9.57433795e-05], f(X) = 0.0000000
[ 97] X = [  4.83771910e-05   8.11690415e-05], f(X) = 0.0000000
[ 98] X = [  3.85394719e-05   6.80521374e-05], f(X) = 0.0000000
[ 99] X = [  2.96855247e-05   5.62469236e-05], f(X) = 0.0000000

array([  2.96855247e-05,   5.62469236e-05])
学習率が低すぎた場合

降下する量が少なすぎて,Xの値がほとんど変わりません

f = lambda X: X[0]**2 + X[1]**2
X = np.array([3.0, 4.0])
gradient_descent(f, X, learning_rate=0.001, max_iter=100)
[  0] X = [ 2.99699995  3.99599995], f(X) = 24.9500243
[  1] X = [ 2.9940029  3.9920039], f(X) = 24.9001485
[  2] X = [ 2.99100885  3.98801185], f(X) = 24.8503724
[  3] X = [ 2.98801779  3.98402378], f(X) = 24.8006958
[  4] X = [ 2.98502972  3.98003971], f(X) = 24.7511185
[  5] X = [ 2.98204464  3.97605962], f(X) = 24.7016403
[  6] X = [ 2.97906255  3.97208351], f(X) = 24.6522611
[  7] X = [ 2.97608343  3.96811138], f(X) = 24.6029805
...
[ 50] X = [ 2.85076078  3.8010152 ], f(X) = 22.5745536
...
[ 96] X = [ 2.72253126  3.63004322], f(X) = 20.5893902
[ 97] X = [ 2.71980868  3.62641313], f(X) = 20.5482314
[ 98] X = [ 2.71708882  3.62278666], f(X) = 20.5071549
[ 99] X = [ 2.71437168  3.61916383], f(X) = 20.4661604

array([ 2.71437168,  3.61916383])
学習率が高すぎた場合

降下する量が多すぎて,Xの値が発散してしまっています

f = lambda X: X[0]**2 + X[1]**2
X = np.array([3.0, 4.0])
gradient_descent(f, X, learning_rate=5.0, max_iter=100)
[  0] X = [-12.00025 -16.00025], f(X) = 400.0140001
[  1] X = [ 48.00075  64.00075], f(X) = 6400.1680012
[  2] X = [-192.00325001 -256.00325001], f(X) = 102402.9120322
[  3] X = [  768.01275031  1024.0127495 ], f(X) = 1638445.6957638
[  4] X = [-3072.05125096 -4096.05124272], f(X) = 26215134.6715417
[  5] X = [ 12288.20484725  16384.20472615], f(X) = 419442142.8762379
[  6] X = [-49152.81880054 -65536.81989439], f(X) = 6711074357.9075851
[  7] X = [ 196611.24598828  262147.29646914], f(X) = 107377187095.1434631
...
[ 50] X = [  1.23730419e+13   1.24100600e+13], f(X) = 307101753022579397827231744.0000000
...
[ 96] X = [  1.23730419e+13   1.24100600e+13], f(X) = 307101753022579397827231744.0000000
[ 97] X = [  1.23730419e+13   1.24100600e+13], f(X) = 307101753022579397827231744.0000000
[ 98] X = [  1.23730419e+13   1.24100600e+13], f(X) = 307101753022579397827231744.0000000
[ 99] X = [  1.23730419e+13   1.24100600e+13], f(X) = 307101753022579397827231744.0000000

array([  1.23730419e+13,   1.24100600e+13])

まとめ

最急降下法アルゴリズムとその実装について紹介しました.

このアルゴリズムは,ニューラルネットワークのパラメタの最適化に利用されています.
具体的には,ニューラルネットワークによる出力を理想状態に近づけるために,実際の出力と理想の出力の誤差を算出する誤差関数の値を最小化するために利用されています.