Juliaで勾配降下法をやってみた

背景

大学の講義でpythonによる機械学習の方法について学んでいるのだが, C言語からプログラミングをはじめて学んだ筆者は, 正直なところpythonのことが嫌いであるそんなに得意ではない.
しかしながら, そんな気持ちを抱えたまま講義を受けるのも面白くないので, 「pythonで受ける講義を, juliaに翻訳してみよう」とふと思い, 挑戦してみることにした.

今回扱う内容：勾配降下法

そもそも, 勾配降下法って？

機械学習を学ぶ上で, 結構最初の方に学ぶであろう 勾配降下法(Gradient Descent)について今回は扱う.
勾配降下法について, IBMのHPでは次のように説明されている.

機械学習モデルやニューラル・ネットワークのトレーニングによく使用される最適化アルゴリズムで, 予測結果と実際の結果の間の誤差を最小限に抑えて機械学習モデルをトレーニングする手法の1つである.

予測結果と実際の結果の誤差はなるべく小さいに越したことは無い. しかしながら, 誤差が小さくなるようなパラメータをやみくもに探してはキリがないので, ある程度山をはって探す必要がある.
勾配降下法はその際に使用される手法の1つなのである.

仕組み

では, どうやって山をはるのか. 勾配降下法ではその名の通り, 勾配を利用して誤差が小さくなるようなパラメータを探すのである.

下の図を見てほしい. 図の横軸$x$はパラメータの値, 縦軸は予測値との誤差を表している.
我々の目的は, 誤差を小さくすることである. この目的を達成するためには, 縦軸が最小となるような$x$, つまり極小値①を取るときの$x$の値を求めることができれば良いのである.

勾配降下法では, まず最初に適当な$x$の値を設定する. (下図のグラフでは仮に$x = 0$としている.)
そうしたら, その地点$x = 0$における傾き, すなわち, $f'(x = 0)$を求める.
この時,

• その地点における傾きが右肩上がり ($f'(x = 0)$が正) であれば求めたい最小値はおそらくもっと左側(つまり, より$x$の値が小さい側)
• 逆に右肩下がり右肩下がり ($f'(x = 0)$が負) であれば求めたい最小値はおそらくもっと右側(つまり, より$x$の値が大きい側)

にあるはずである.

また, 勾配が急であればあるほど, 求めたい$x$の値より遠いはずである.
よって, これらの情報を加味して, $x$の値を次のように更新する.
$$x_{t+1} = x_{t} - \alpha\cdot f'(x=x_{t})
$$
ここで, $\alpha$は学習率といい, この値が大きければ大きいほど, 値がより大きく更新されることになる.
今回の例では, $x = 0$を仮の最初の値としたため, 0を起点に$x$のパラメータを更新する.

そして更新された値$x'$をもとに, 再度「勾配を調べる→パラメータを更新→…」という作業を繰り返し, 最終的に求めたいパラメータの値にたどり着く..という訳である.

ここで注意点が2つある.
1点目は, 適切な学習率$\alpha$を設定しないと, うまく収束しない場合がある点である.
前述の通り, 学習率が大きければ大きいほど, 一回の作業でパラメータがより大きく更新されることになる.
ここで誤ってパラメータを大きくし過ぎると, 極小値を取るパラメータを大きく通りすぎてしまう危険性がある.　(下図の左側のイメージ)
2点目は, 初期値によって誤った極小値が求まってしまう場合がある点である.
下図の右側の例がまさしくそのパターンであるが, 勾配降下法はその時点における傾きをもとにパラメータを更新するため, もしも, 初期位置の近傍に局所的な極小値②が存在していた場合, そちらの方に吸い込まれてしまう危険性がある.

これを回避する方法は他にもあると思うが, 筆者にはまだそれを言及できるほど知識がないため, ここでは省略させていただく. (いずれは書きたい.)

ではこれらの知識をもとに, Juliaでの実装を行ってみようと思う.

Juliaによる実装

今回は$f(x) = x^2 - 10x - 80$の極小値を求めてみようと思う.
まずは, お題の関数の定義である.

# お題となる関数
function my_function(x)
    x^2 - 10x - 80
end

続いて, 微分の機能を果たす関数の定義である.

# 数値微分を定義する関数
function numeric_gradient(f, x, h=1e-5)
    (f(x + h) - f(x - h)) / (2 * h)
end

ここで, プログラミングにおいては真の意味での微分(解析的微分)を行うことができないため, 微小区間h=1e-5における傾きを調べる方法で対処するとする. このように具体的な数値を用いて微分を行う手法を数値的微分という.

そして次に, 本題である勾配降下法の定義である.

# 勾配降下法を用いた最適化
function gradient_descent(initial_x, learning_rate, num_iterations)
    x = initial_x
    x_history = Array{Float64}(undef, num_iterations)
    x_history[begin] = x
    for i in 2:num_iterations
        gradient = numeric_gradient(my_function, x)
        x = x - gradient * learning_rate
        x_history[i] = x
        if i % 10 == 0
            println("x=$(x_history[i])のとき, y=$(my_function(x_history[i]))" )   
        end
    end
    return x_history
end

initial_xは$x$の初期値, learning_rateは学習率$\alpha$, num_iterationsはパラメータの更新回数を表す.
x_historyは, 更新されてゆくパラメータ$x$の履歴を記録しておくもので, リストとなっている. (JuliaはFortranのように配列のインデックスが1から始まる点に注意.)

ここで,

for i in 2:num_iterations
    gradient = numeric_gradient(my_function, x)
    x = x - gradient * learning_rate
    x_history[i] = x
    # 10回繰り返すごとに, f(x)の値を求める
    if i % 10 == 0
        println("x=$(x_history[i])のとき, y=$(my_function(x_history[i]))" )   
    end
end

に注目してみる.

gradient = numeric_gradient(my_function, x)

はまさしく, $x$における傾きをnumeric_gradient関数によって求め, 変数gradientに代入している.
また, その次の行の

x = x - gradient * learning_rate

では, 求めた傾きと学習率$\alpha$を用いて, パラメータを更新している.

これらをつなげると, 次のようなコードになる.

# 極小値を求める関数
function my_function(x)
    x^2 - 10x - 80
end

# 数値微分を定義する関数
function numeric_gradient(f, x, h=1e-5)
    (f(x + h) - f(x - h)) / (2 * h)
end

# 勾配降下法を用いた最適化
function gradient_descent(initial_x, learning_rate, num_iterations)
    x = initial_x
    x_history = Array{Float64}(undef, num_iterations)
    x_history[begin] = x
    for i in 2:num_iterations
        gradient = numeric_gradient(my_function, x)
        x = x - gradient * learning_rate
        x_history[i] = x
        # 10回繰り返すごとに, f(x)の値を求める
        if i % 10 == 0
            println("x=$(x_history[i])のとき, y=$(my_function(x_history[i]))" )   
        end
    end
    return x_history
end

# 実際に求める
initial_x = 0 # 初期xの定義
learning_rate = 0.05 # 学習率
num_iterations = 100 # 繰り返し回数
x_history = gradient_descent(initial_x, learning_rate, num_iterations)

これの実行結果は以下のとおりである.

x=3.062897555032862のとき, y=-101.24763411770233
x=4.324574141243431のとき, y=-104.54379990932296
x=4.76449356511921のとき, y=-104.94453671912974
x=4.917883983672766のとき, y=-104.99325695986255
x=4.971367915516112のとき, y=-104.9991802037381
x=4.990016609518477のとき, y=-104.99990033191449
x=4.996519006965627のとき, y=-104.99998788268749
x=4.99878625277006のとき, y=-104.99999852681766
x=4.999576792492633のとき, y=-104.9999998208954
x=4.999852436675667のとき, y=-104.99999997822506

実行結果を見ると, おおよそ$x=5$で収束していることが分かる. 今回極小値を求める関数も
$$
\begin{aligned}
f(x) &= x^2 - 10x - 80 \\
&= (x-5)^2 -105
\end{aligned}
$$
なので, 頂点が$(5, -105)$となる下に凸の二次関数グラフであることより, うまく極小値が求められていることが分かるだろう.

下記のコードを用いて, グラフに直してみると,

# 極小値を求めていることの可視化
using PyPlot

fig, ax = subplots(figsize=(8,3))

# 関数の描写をするための値
xs = range(0, 6, length=100)

ax.plot(xs, my_function.(xs), "b")
ax.plot(x_history, my_function.(x_history), "rx", label="sequence")

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.grid()
ax.legend()

display(fig)

となり, 実際に極小値を求めることができていることが分かる.