Pynote

Python、機械学習、画像処理について

数学 - 勾配法について可視化して理解する。

最適化

与えられた制約条件のもとで関数の値を最大化または最小化する変数の値を求めることを最適化という。

勾配法

勾配法の仕組み

\mathbb{R}^n の開集合 \Omega 上で n 次元関数 f: \Omega \to \mathbb{R} が定義されているとする。
この関数の極小値を求めることを考える。

今、ある点 \boldsymbol{x}_k にいる場合、その点の勾配 \nabla f (\boldsymbol{x}_k) は点 \boldsymbol{x}_k の近傍で関数の値を最も増加させる「最急上昇方向 (steepest ascent)」であるから、その反対方向 -\nabla f (\boldsymbol{x}_k) は関数の値を最も減少させる「最急降下方向 (steepest descent)」である。(証明は以下の記事を参照)

pynote.hatenablog.com

よって、点 \boldsymbol{x}_k から勾配と逆方向に進んで、新しい点 \boldsymbol{x}_{k + 1} に移動することで関数の値を減少させられる。
 \displaystyle
\boldsymbol{x}_{k + 1} = \boldsymbol{x}_k -\alpha \nabla f (\boldsymbol{x}_k)

上記の更新を1回行うことを「1反復 (iteration)」という。

\alpha \in \mathbb{R}^+ は勾配方向にどのくらい進むかを表すスカラーであり、ステップ幅という。
各反復で関数値が減少するよう移動したいので、f(\boldsymbol{x}_k) \ge f(\boldsymbol{x}_{k + 1}) を満たす \alpha を選択する。
ステップ幅は反復ごとに変更してよい。

反復は点が極小点に到達するまで繰り返す。

極小点に到達したかどうかは、極値点では勾配が \nabla f (\boldsymbol{x}_k) = \boldsymbol{0} となることを利用して判定する。
数値計算においては、厳密に \boldsymbol{0} になることはないので、十分 \boldsymbol{0} に近ければ収束したと判定し、反復を終了する。

これまでの一連の流れを整理すると、以下のようになる。

[アルゴリズム] 最急降下法

1. 初期点 \boldsymbol{x}_0 を決める。
2. 勾配 \nabla f (\boldsymbol{x}_k) を計算し、\nabla f (\boldsymbol{x}_k) = \boldsymbol{0} の場合、終了する。
3. 点を次のように更新する。
 \displaystyle
\boldsymbol{x}_{k + 1} = \boldsymbol{x}_k -\alpha \nabla f (\boldsymbol{x}_k)
4. 2 に戻る。

また、n 次元関数 f: \Omega \to \mathbb{R} の極大点を求める勾配法は「最急上昇法」と言われる。

[アルゴリズム] 最急上昇法

1. 初期点 \boldsymbol{x}_0 を決める。
2. 勾配 \nabla f (\boldsymbol{x}_k) を計算し、\nabla f (\boldsymbol{x}_k) = \boldsymbol{0} の場合、終了する。
3. 点を次のように更新する。
 \displaystyle
\boldsymbol{x}_{k + 1} = \boldsymbol{x}_k + \alpha \nabla f (\boldsymbol{x}_k)
4. 2 に戻る。

ステップ幅の決め方

勾配法ではステップ幅 \alpha をどのように決めるかが問題になる。
値が小さすぎると収束までに時間がかかってしまい、逆に大きすぎるといつまでの収束しないか発散してしまう。

ステップ幅を直線探索で決める。

1つの方法として、探索方向で関数値が最小となるようにステップ幅を選ぶ方法が考えられる。これを直線探索 (line search)という。
今、点 \boldsymbol{x} にいる場合、探索方向の直線上の点は \boldsymbol{x}'(\alpha) = \boldsymbol{x} -\alpha \nabla f (\boldsymbol{x}), \  \alpha \in \mathbb{R}^+ と表せる。これを探索直線 (search line)という。
数式で表すと、

 \displaystyle
\min_{\alpha} f(\boldsymbol{x}'(\alpha))

を満たす \alpha を探せばよい。
これを満たす点は極小点なので、F(\alpha) = f(\boldsymbol{x}'(\alpha)) とおくと、\frac{dF(\alpha)}{d\alpha} = 0 を満たす。

[定理] 直線探索でステップ幅を決めた場合、勾配は1つ前の勾配と直交する。

\boldsymbol{x} = (x_1, x_2, \cdots, x_n)^T とすると、関数 Fx_1, x_2, \cdots, x_n によって決まり、x_1, x_2, \cdots, x_n\alpha によって決まるので連鎖律が適用でき、

 \displaystyle
\frac{dF(\alpha)}{d\alpha} = \sum_{i=1}^{n} \frac{\partial f}{\partial x'_i} \frac{\partial x'_i}{\partial \alpha}

x'_i(t) = x_i + \alpha \frac{\partial f}{\partial x} であるから、\frac{\partial x'_i}{\partial \alpha} = \frac{\partial f}{\partial x} となる。
よって、先程の式に代入すると、

 \displaystyle
\frac{dF(t)}{dt} = \sum_{i=1}^{n} \frac{\partial f}{\partial x'} \frac{\partial f}{\partial x_i}
= (\nabla f(\boldsymbol{x}'), \nabla f(\boldsymbol{x}))

よって、探索直線上で極小値をとる点では、(\nabla f(\boldsymbol{x}), \nabla f(\boldsymbol{x}')) = 0 を満たす。
つまり、探索方向 \nabla f(\boldsymbol{x}) と探索直線の極小値をとる点での勾配 \nabla f(\boldsymbol{x}') が直交することを意味している。

探索直線の可視化

2変数関数 f(x_1, x_2) = x_1^2 + x_2^2 + x_1 x_2 を定義する。
この勾配は \nabla f(x_1, x_2) = \left(2 x_1 + x_2, 2 x_2 + x_1\right)^T である。

関数を可視化する。

関数 f を描画する。

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import Axes3D, proj3d


def f(x):
    x1, x2 = x
    return x1**2 + x2**2 + x1 * x2


def f_prime(x):
    x1, x2 = x
    return np.array([2 * x1 + x2, 2 * x2 + x1])


X1, X2 = np.mgrid[-10:11, -10:11]
Y = f((X1, X2))  # 各点での関数 f の値を計算する。

fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(111, projection='3d')

# 各軸のラベルを設定する。
ax.set_xlabel('$x$', fontsize=15)
ax.set_ylabel('$y$', fontsize=15)
ax.set_zlabel('$z$', fontsize=15)

# グラフを作成する。
ax.plot_surface(X1, X2, Y, alpha=0.4, antialiased=False)

# 視点を設定する。
ax.view_init(elev=50, azim=120)


勾配および探索方向を可視化する。

今、点 (7, 3) にいるとする。
(7, 3) 及び対応する関数値 f(7, 3) を描画する。

x = 7, 3
y = f(x)
delta = f_prime(x)

# 点 (3, 7) 及び関数値 f(3, 7) を作成する。
ax.scatter(*x, 0, color='DarkBlue', s=50)  # 点 (3, 7)
ax.scatter(*x, y, color='DarkSlateGray', s=50)  # f(3, 7)
ax.text(x[0], x[1] - 1, y, 'f({}, {})'.format(*x), size=16)  # 点 (3, 7) のラベル
ax.text(x[0], x[1] + 4, 0, '({}, {})'.format(*x), size=16)  # f(3, 7) のラベル
# 点 (3, 7) 及び関数値 f(3, 7) を結ぶ線
ax.plot([x[0], x[0]], [x[1], x[1]], [0, y], c='black', lw=2)

(7, 3) における勾配は \nabla f(7, 3) = (17, 13)^T である。
最急上昇方向 \nabla f(7, 3) = (17, 13)^T、最急降下方向 -\nabla f(7, 3) = (17, 13)^T を描画する。

class Arrow3D(FancyArrowPatch):
    def __init__(self, xs, ys, zs, *args, **kwargs):
        FancyArrowPatch.__init__(self, (0, 0), (0, 0), *args, **kwargs)
        self._verts3d = xs, ys, zs

    def draw(self, renderer):
        xs3d, ys3d, zs3d = self._verts3d
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
        self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
        FancyArrowPatch.draw(self, renderer)


# 最急降下方向
descent_head = x + delta * 0.2
ax.text(descent_head[0], descent_head[1], 0, 'descent',
        size=16, horizontalalignment='right')
descent = Arrow3D([x[0], descent_head[0]], [x[1], descent_head[1]], [0, 0], mutation_scale=10,
                  lw=3, arrowstyle='-|>', color='r')
ax.add_artist(descent)

# 最急上昇方向
ascent_head = x - delta * 0.2
ax.text(ascent_head[0], ascent_head[1] + 2, 0,
        'ascent', size=16, horizontalalignment='left')
ascent = Arrow3D([x[0], ascent_head[0] + 0.5], [x[1], ascent_head[1]], [0, 0], mutation_scale=10,
                 lw=3, arrowstyle='-|>', color='r')
ax.add_artist(ascent)


探索直線を可視化する。

探索直線上の点は \nabla (3, 7) - \alpha (17, 13)^T と表せる。
青の直線が探索直線、緑の線が探索直線上の各点に対応する関数値を表している。

# [-10, 10]^2 の範囲の直線上の点
line_X = np.array([x - step * delta for step in np.arange(-1, 1, 0.01)])
line_X = line_X[np.all((-10 <= line_X) & (line_X <= 10), axis=1)]
# 直線上の点に対応する関数値
line_Y = f(line_X.T)

# 直線及び対応する関数値を作成する。
ax.plot(line_X[:, 0], line_X[:, 1], line_Y, c='g')
ax.plot(line_X[:, 0], line_X[:, 1], 0, c='b')


直線探索でステップ幅を決める場合の勾配法

勾配法の結果を可視化するための関数を定義しておく。

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.lines import Line2D
from mpl_toolkits.mplot3d import Axes3D


def draw_history(f, xs, ys, elev=70, azim=-70):
    fig = plt.figure(figsize=(15, 7))

    X1, X2 = np.mgrid[-10:11, -10:11]
    Y = f((X1, X2))  # 各点での関数 f の値を計算する。

    # 勾配のベクトル図を作成する。
    ax1 = fig.add_subplot(121)
    ax1.set_title('Contour')
    ax1.set_xlabel('$x$', fontsize=15)
    ax1.set_ylabel('$y$', fontsize=15)
    ax1.plot(xs[:, 0], xs[:, 1], 'ro-', mec='b', mfc='b', ms=4)
    contours = ax1.contour(X1, X2, Y)
    ax1.clabel(contours, inline=1, fontsize=10, fmt='%.2f')

    # グラフを作成する。
    ax2 = fig.add_subplot(122, projection='3d')
    ax2.set_title('Surface')
    ax2.set_xlabel('$x$', fontsize=15)
    ax2.set_ylabel('$y$', fontsize=15)
    ax2.set_zlabel('$z$', fontsize=15)
    ax2.plot(xs[:, 0], xs[:, 1], ys,
             'ro-', mec='b', mfc='b', ms=4)
    ax2.plot_surface(X1, X2, Y, alpha=0.3, edgecolor='black')
    ax2.view_init(elev=elev, azim=azim)

    plt.show()

勾配法を行う関数を実装する。
数値計算丸め誤差があるため、勾配の値が厳密に \boldsymbol{0} にならない。
そのため、勾配の各要素と \boldsymbol{0} の差が 0.001 未満であれば、収束したと判断している。
直線探索を行い、ステップ幅を決めるには scipy.optimize.line_search() を使う。

from scipy.optimize import line_search

def gradient_decent(f, f_prime, init_x, max_iters=100):
    x = np.array(init_x, dtype=float)  # 初期点
    xs, ys = [], []  # x, f(x) の履歴

    for itr in range(1, max_iters + 1):
        xs.append(x), ys.append(f(x))

        delta = f_prime(x)  # 勾配を計算する。
        if np.all(np.abs(delta) < 0.001):
            break  # 極値に到達

        # 直線探索でステップ幅を決定する。
        step = line_search(f, f_prime, xk=x, pk=-delta)[0]
        x = x - step * delta  # 更新する。

    return itr, np.array(xs), np.array(ys)
# 関数及び勾配を定義する。
def f(x):
    x1, x2 = x
    return x1**2 + x2**2 + x1 * x2


def f_prime(x):
    x1, x2 = x
    return np.array([2 * x1 + x2, 2 * x2 + x1])


# 勾配法を実行する。
num_itrs, xs, ys = gradient_decent(f, f_prime, init_x=np.array([7., 3.]))
print('iterations', num_itrs)  # iterations 4

# 結果を可視化する。
draw_history(f, xs, ys)

結果を見ると、先程の定理の通り、ある反復の勾配とその前後の反復の勾配が直交していることが確認できる。

ステップ幅を定数で決める。

反復ごとに直線探索を行うのは、計算量が多くなってしまう。
その代わりにステップ幅を適当な小さい値に設定する方法が一般的に用いられている。
この定数の適切な値は関数の性質によって変わってくる。
以下、ステップ幅を定数 0.1 とした場合の勾配法を実装する。

def gradient_decent(f, f_prime, init_x, step=0.1, max_iters=100):
    x = np.array(init_x, dtype=float)  # 初期点
    xs, ys = [], []  # x, f(x) の履歴

    for itr in range(1, max_iters + 1):
        xs.append(x), ys.append(f(x))

        delta = f_prime(x)  # 勾配を計算する。
        if np.all(np.abs(delta) < 0.001):
            break  # 極値に到達

        x = x - step * delta  # 更新する。

    return itr, np.array(xs), np.array(ys)
# 関数及び勾配を定義する。
def f(x):
    x1, x2 = x
    return x1**2 + x2**2 + x1 * x2


def f_prime(x):
    x1, x2 = x
    return np.array([2 * x1 + x2, 2 * x2 + x1])


# 勾配法を実行する。
num_itrs, xs, ys = gradient_decent(f, f_prime, init_x=np.array([7., 3.]))
print('iterations', num_itrs)  # iterations 74

# 結果を可視化する。
draw_history(f, xs, ys)


勾配法の欠点について

関数の形状によっては収束に必要な反復回数が増えてしまう。

例えば、2変数関数 f(x_1, x_2) = 0.05 x_1^2 + x_2^2 は以下の画像のような谷の形をしている。
このような関数に対して、勾配法を実行すると、谷と垂直の方向の勾配はなだらかのため、極値点に到達するまでに多くの反復回数が必要となる。

# 関数及び勾配を定義する。
def f(x):
    x1, x2 = x
    return 0.05 * x1**2 + x2**2


def f_prime(x):
    x1, x2 = x
    return np.array([0.1 * x1, 2 * x2])


# 勾配法を実行する。
num_itrs, xs, ys = gradient_decent(
    f, f_prime, step=0.9, init_x=(7, 3))
print('iterations', num_itrs)  # iterations 71

# 結果を可視化する。
draw_history(f, xs, ys)


局所解に収束してしまう。

収束判定条件が勾配が \boldsymbol{0} の場合であるため、極値が複数存在する場合、最小値でない極小値に収束しまう場合がある。
例えば、2変数関数 f(x_1, x_2) = \sin(x_1 + x_2) + 0.5 (x_1 + x_2) は以下の波状の関数である。
これは極小値が複数存在するため、勾配法を実行した際にそこにハマってしまいっている。

# 関数及び勾配を定義する。
def f(x):
    x1, x2 = x
    return np.sin(x1 + x2) + 0.5 * (x1 + x2) + 1  # McCormick function


def f_prime(x):
    x1, x2 = x
    return np.array([np.cos(x1 + x2) + 0.5, np.cos(x1 + x2) + 0.5])


# 勾配法を実行する。
num_itrs, xs, ys = gradient_decent(
    f, f_prime, step=0.5, init_x=(9, 4))
print('iterations', num_itrs)

# 結果を可視化する。
draw_history(f, xs, ys, 20, -50)


鞍点に収束してしまう。

収束判定条件が勾配が \boldsymbol{0} の場合であるため、鞍点または停留点に収束してしまう場合がある。
例えば、2変数関数 f(x_1, x_2) = x_1^2 + x_2^2 + x_1 x_2^2 (0, 0)^T の鞍点を持つ関数である。
鞍点で勾配の値が \boldsymbol{0} となるため、勾配法を実行した際にそこにハマってしまいっている。

# 関数及び勾配を定義する。
def f(x):
    x1, x2 = x
    return x1**2 + x2**2 + x1 * x2 ** 2


def f_prime(x):
    x1, x2 = x
    return np.array([2 * x1 + x2, 2 * x2 + x1])


# 勾配法を実行する。
num_itrs, xs, ys = gradient_decent(
    f, f_prime, step=0.5, init_x=(9, 4))
print('iterations', num_itrs)

# 結果を可視化する。
draw_history(f, xs, ys)