最急降下法（勾配法）のアルゴリズム
ある関数の最小値を数値的に求める方法として、古くから最急降下法（勾配降下法）というアルゴリズムが知られています。最急降下法は適当な所にボールを置いて手を離し、斜面に沿って転がしながら、最終的に落ち着いた場所を最小値とするようなイメージのアルゴリズムです。
簡単な例として、二次関数 $f(x)=x^2$ の最小値を最急降下法で求めてみましょう。
この関数の点 $x$ における勾配（傾き）は導関数 $f'(x)=2x$ で与えられます。
最初に適当な点（初期値）$x_0$ を決めます。
どこでもいいのですが、とりあえず $x_0=2$ としておきます。
初期値における勾配は
\[f'(x_0)=2x_0=4\tag{1}\]
です。上って行く方向とは逆に点を動かすために、この勾配に学習率とよばれる正の係数 $\alpha$ を掛けた値を $x_0$ から引いて（すなわち $x$ 軸の負の方向に点を移動させて）、その点を $x_1$ とします。
\[x_1=x_0-\alpha f'(x_0)=2-4\alpha\tag{2}\]
$\alpha=0.1$ と決めると
\[x_1=2-0.4=1.6\tag{3}\]
となって、少しだけ原点に近づいたことがわかります。この点で再び勾配を計算します。
\[f'(x_1)=2\times 1.6=3.2\tag{4}\]
先ほどと同じように、この勾配を使って次の点 $x_2$ を決めます。
\[x_2=x_1-\alpha f'(x_1)=1.6-0.1\times 3.2=1.6-0.32=1.28\tag{5}\]
また少し原点に近づきました。さらに、これを繰り返して、$x_3,\ x_4,\ ……$ を求めていけば、いずれは最小値をとる場所（原点）に近い所に達します。一般化すると、$x_k$ は
\[x_k=x_{k-1}-\alpha f'(x_{k-1})\tag{6}\]
によって与えられます。以下のコードでは上記のアルゴリズムにしたがって最小値を探索しながら、Matplotlib を使って初期値から窪みの底まで移動する点の軌跡をプロットします。
# Python Steepest descent
# In[1]
# 最急降下法による最小値探索[1]
# NumPyをインポート
import numpy as np
# matplotlib.pyplotをインポート
import matplotlib.pyplot as plt
# 配列要素を小数点以下8桁まで表示
np.set_printoptions(precision=8, floatmode="fixed")
# 2次関数の定義
def func(x):
return x**2
# func(x)の導関数と導関数の絶対値
def dfunc(x):
df = 2*x
df_abs = np.abs(df)
return df, df_abs
# ★★★ 2次関数グラフの描画 ★★★
# Figureを作成してAxesを１つ追加
fig = plt.figure()
ax = fig.add_subplot(111)
# Axesのタイトルを設定
ax.set_title("y = x**2", fontsize=16)
# 軸ラベルと軸範囲を設定
ax.set_xlabel("x", fontsize=16)
ax.set_ylabel("y", fontsize=16)
ax.set_xlim([-4, 4])
ax.set_ylim([0, 16])
# -4～4まで0.1刻みの数値の配列
x = np.arange(-4, 4, 0.1)
# Axesに2次関数のグラフを灰色で描画
# 描画優先順位zorder=1
ax.plot(x, func(x), color="gray", zorder=1)
# ★★★ 最急降下法による最小値の探索 ★★★
# xの初期値
x = 2
# 学習率α
alpha = 0.01
# 収束判定条件
eps = 0.001
# 繰返しの最大数
k_max = 1000
for k in range(1, k_max):
x -= alpha * dfunc(x)[0]
# (x, f(x))を散布図にプロット
# マーカーサイズs=5,描画優先順位：zorder=2
ax.scatter(x, func(x), s=5, color="red", zorder=2)
# 勾配の絶対値がeps以下になったらループ終了
if dfunc(x)[1] < eps:
break
# 最小値をとるxと勾配の大きさ
result = np.array([x, dfunc(x)[1]])
print(result)
# [0.00049539 0.00099078]
実行結果を見ると、$x=0$ に近い値が計算されています。
収束判定条件 eps をもっと小さくすれば精度を上げることができます。グラフを見ると、初期値 $x=2$ から $x=0$ まで点が移動していく様子がわかります。
$\alpha$ を大きくとれば、1 step の移動距離を増やすことができますが、あまり大きすぎると、最小値のある点を大きく超えて反復運動を始めてしまい、最小値に落ち着くまでかえって時間がかかることもあります。
複数の極小値をもつ関数の最小値探索
最急降下法は初期値の近くの窪みに落ち込んでしまうので、最小値というよりは、極小値を探すアルゴリズムです。複数の極小値をもつ関数に最急降下法を適用すると、$x$ の初期値の選び方によっては浅い窪みに落ちて抜け出せなくなってしまう可能性があります。このような関数を扱う場合は、複数の初期値を選んで極小値となる $x$ を計算させて $f(x)$ の値を比較します。
例として、$f(x)=x\cos 2x$ という関数の $0 \leq x \leq 2\pi$ の範囲における最小値を見つけます。下のグラフを見て分かる通り、この関数は範囲内に 2 つの極小値をもっています。
$f(x)$ を微分すると、$f'(x)=\cos 2x-2x\sin 2x$ となります。
初期値として $x_0=1,\ 2,\ 3,\ 4,\ 5$ の点を選んで、それぞれの地点から点を転がして極小値を計算してみます。
# In[2]
# 最急降下法による最小値探索[2]
# 関数の定義
def func(x):
return x * np.cos(2*x)
# 導関数
def dfunc(x):
df = np.cos(2*x) - 2*x*np.sin(2*x)
df_abs = np.abs(df)
return df, df_abs
# ★★★ グラフの描画 ★★★
# Figureを作成してAxesを１つ追加
fig = plt.figure()
ax = fig.add_subplot(111)
# Axesのタイトルを設定
ax.set_title("y = xcos2x", fontsize=16)
# 軸ラベルと軸範囲を設定
ax.set_xlabel("x", fontsize=16)
ax.set_ylabel("y", fontsize=16)
ax.set_xlim([0, 2*np.pi])
# 0～2piまで0.1刻みの数値の配列を定義
x = np.arange(0, 2*np.pi, 0.1)
# Axesにグラフをプロット
# 描画優先順位：zorder=2
ax.plot(x, func(x), color="gray", zorder=1)
# ★★★ 最急降下法による最小値の探索 ★★★
# 学習率α
alpha = 0.01
# 収束判定条件
eps = 0.001
# 繰返しの最大数
k_max = 1000
# 初期値=1,2,3,4,5
for j in range(1, 6):
x = j
for k in range(1, k_max):
x -= alpha * dfunc(x)[0]
# (x, f(x))を散布図にプロット
# マーカーサイズ：s=5,描画優先順位：zorder=2
ax.scatter(x, func(x), s=5, color="red", zorder=2)
# 勾配の絶対値がeps以下になったらループ終了
if dfunc(x)[1] < eps:
break
# x,f(x)の値を表示
print(np.array([x, func(x)]))
# [ 1.71268122 -1.64418563]
# [ 1.71293651 -1.64418564]
# [ 1.71293529 -1.64418564]
# [ 4.76462398 -4.73864711]
# [ 4.76471852 -4.73864710]
実行結果を見ると、2 つの極小値が発見されていて、$x=4.765$ あたりで最小値となることがわかります。グラフを見ると、初期値 $x_0=1,\ 2,\ 3$ からスタートした点は浅い窪みに、$x_0=4,\ 5$ からスタートした点は深い窪みに落ちていることがわかります。
二変数関数(曲面)の最小値探索
次は二変数関数 $z=f(x, y)$ の最小値を最急降下法で計算してみましょう。曲面上の各点 $(x, y)$ における勾配 (gradient) は
\[\nabla f=\begin{pmatrix}\cfrac{\partial f}{\partial x}\\\cfrac{\partial f}{\partial y}\end{pmatrix}\tag{7}\]
というベクトルで表され、その大きさは
\[|\nabla f|=\sqrt{\left(\frac{\partial f}{\partial x}\right)^2+\left(\frac{\partial f}{\partial y}\right)^2}\tag{8}\]
で計算できます。したがって、あるステップにおける点の位置は学習率 $\alpha$ を用いて
\[\begin{pmatrix}x_k\\ y_k\end{pmatrix}=\begin{pmatrix}x_{k-1}\\ y_{k-1}\end{pmatrix}-\alpha\begin{pmatrix}\left.\cfrac{\partial f}{\partial x}\right|_{x_{k-1}}\\\left.\cfrac{\partial f}{\partial y}\right|_{x_{k-1}}\end{pmatrix}\tag{9}\]
と表すことができます。例として $z=x^2+y^2$ の最小値を求めるコードを書いてみます。この関数の勾配を求めておきましょう。$z$ を $x$ と $y$ で偏微分すると
\[\frac{\partial z}{\partial x}=2x,\quad \frac{\partial z}{\partial y}=2y\tag{10}\]
となるので、この関数の各点における勾配は $\displaystyle\binom{2x}{2y}$、その大きさは $\sqrt{(2x)^2+(2y)^2}$ で与えられます。
# In[3]
# 最急降下法による最小値探索[3]
# 関数の定義
def func(x ,y):
return x**2 + y**2
# 勾配と勾配の大きさ(ノルム)
def grad(x, y):
gf = np.array([2*x, 2*y])
gf_norm = np.linalg.norm(gf)
return gf, gf_norm
# Figureを追加
fig = plt.figure(figsize = (10, 6))
# 3DAxeSを追加
ax = fig.add_subplot(111, projection="3d")
# 軸ラベルを設定
ax.set_xlabel("x", size=16)
ax.set_ylabel("y", size=16)
ax.set_zlabel("z", size=16)
# 格子点の作成
X, Y = np.mgrid[-10:10, -10:10]
# 高度の計算式
Z = func(X, Y)
# ワイヤーフレームを描画
ax.plot_wireframe(X, Y, Z, color="blue", zorder=1)
# 視点の設定
ax.view_init(45, 45)
# ★★★ 最急降下法による最小値の探索 ★★★
# 初期値の設定
x, y = (-8, -4)
# 学習率α
alpha = 0.01
# 収束判定条件
eps = 0.001
# 繰返しの最大数
k_max = 10000
for k in range(1, k_max):
x -= alpha * grad(x, y)[0][0]
y -= alpha * grad(x, y)[0][1]
# (x, f(x))を散布図にプロット
ax.scatter(x, y, func(x, y), s=2, color="red", zorder=2)
# 勾配の大きさがeps以下になったらループ終了
if grad(x, y)[1] < eps:
break
# 最小値をとるx,yを表示
print(np.array([x, y]))
# [-0.00044436 -0.00022218]
曲面はワイヤーフレーム、点の軌跡は散布図で描いています。
初期値から窪みに落ちてゆく様子がわかります。
最終的には、ほぼ $(0, 0)$ の点に落ち着いています。
細かい事を言えばきりがないのですが、中には勾配がどの点においても一定であるような関数も存在していて（たとえば $z=\sqrt{x^2+y^2}$）、そのような関数では窪みが尖っているので、極小値の近くに達しても勾配は一向に緩やかにならず、eps の判定にかかりません。それでも学習率の値を小さくして、繰り返し回数を十分に大きくとっておけば、最急降下法のアルゴリズムにしたがって点は下へ下へと落ちていくので、最小値を発見することはできます。多少はステップごとに値が揺れますが、学習率の値が小さければ、窪みの中心付近における反復運動はとても小さくなります。
以下の記事は補足的内容なので、学習をお急ぎの方は【次の記事】へそのままお進みください。
転移機能付最小値探索関数
任意の関数と偏導関数を与えて最小値を探索する steepest() 関数のコードを載せておきます。この関数は最小値以外の極小値への落ち込みをランダム転移によって避けるように工夫されています（ほぼ確実に最小値を見つけ出します）。少し長いコードですが、わからない箇所があれば遠慮なくコメント欄で質問をお寄せください。
# In[4]
# 最急降下法による最小値探索関数
def steepest(axes, func, grad, xrange, yrange,
alpha=0.01, eps=0.001, k_max=1000, plot=True):
# 格子点の作成
X, Y = np.mgrid[xrange[0]:xrange[1], yrange[0]:yrange[1]]
# 高度の計算式
Z = func(X, Y)
# 等高線を描画
if plot == True:
# 等高線をプロット
ct = axes.contour(X, Y, Z, 12, colors="black")
axes.clabel(ct, fmt="%0.1f", fontsize=10)
# x,yの初期値
x = np.random.uniform(xrange[0], xrange[1])
y = np.random.uniform(yrange[0], yrange[1])
for k in range(1, k_max):
x -= alpha * grad(x, y)[0][0]
y -= alpha * grad(x, y)[0][1]
# 10ステップごとに無作為に選んだ点と現在の点の高度を比較して
# 無作為に選んだ点のほうが低ければ転移する
if k % 10 == 0:
x2 = np.random.uniform(xrange[0], xrange[1])
y2 = np.random.uniform(yrange[0], yrange[1])
if func(x2, y2) < func(x, y):
x = x2
y = y2
# 指定範囲を超えそうになったら位置をリセット
if not (xrange[0] < x < xrange[1]
and yrange[0] < y < yrange[1]) :
x = np.random.uniform(xrange[0], xrange[1])
y = np.random.uniform(yrange[0], yrange[1])
# 軌跡を等高線にプロット
if plot == True:
axes.plot(x, y, color="red",
marker="o", markersize=1,
linestyle="none")
# 勾配の大きさがeps以下になったらループ終了
if grad(x, y)[1] < eps:
break
return np.array([x, y])
steepest() は関数の凸凹がわかりやすいように、曲面の代わりに等高線が表示されるようになっています。試しに二変数関数
\[z(x,y)=\sin x+\sin\sqrt{x}+\sin\frac{y}{2}\tag{11}\]
の最小値を見つけてみましょう。この関数の $x,\ y$ による偏微分は
\[z_x=\cos x+\frac{\cos\sqrt{x}}{2\sqrt{x}}\quad z_y=\frac{1}{2}\cos y\tag{12}\]
となっています。
# In[5]
# 最急降下法による最小値探索[4]
# 疑似乱数シードを設定
np.random.seed(11)
# 関数の定義
def func2(x ,y):
return np.sin(x) + np.sin(np.sqrt(x)) + np.sin(y/2)
# 勾配と勾配のノルム
def grad2(x, y):
fx = np.cos(x) + 0.5 * np.cos(np.sqrt(x))/(np.sqrt(x))
fy = 0.25 * np.cos(y/2)
gf = np.array([fx, fy])
gf_norm = np.linalg.norm(gf)
return gf, gf_norm
# FigureとAxesの設定
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111)
ax.set_xlabel("x", size=16)
ax.set_ylabel("y", size=16)
# 探索範囲
xr = [0, 15]
yr = [0, 15]
# 最急降下法関数で指定範囲内の最小値を探索
s = steepest(ax, func2, grad2, xr, yr, alpha=0.2)
print(s)
# [11.14300987 9.41684212]
指定範囲内に浅い窪みと深い窪み (最小値) があります。最初は浅いくぼみのほうへ落ちそうになりますが、ランダム転移機能によって抜け出して、最終的には深い窪みの底に達していることがわかります。等高線が表示されたほうがわかりやすいと思いますが、それなりに処理時間を要します。plot 引数を False に設定すると、軌跡のプロットは行われずに最小値をとる x と y だけを返します。
コメント
下記は誤植と思われますので、ご確認ください。
式番号で、(8) が 2 ヶある。
In[4]プログラムで、# ワイヤフレームを描画 → # 等高線を描画
In[5]プログラムの下の文章で、最小値だけ → 最小値をとる x, y だけ
ありがとうございます。
直しておきました。