Python
機械学習
数学
最適化
人工知能

機械学習のためのニュートン法 (1 変数から多変数まで)

はじめに

本記事では,金谷健一先生の「これなら分かる最適化数学1 を参考に,ニュートン法を題材にして少し数学的な復習をしつつ,Python を使って最適化ってこんな感じか〜というのを眺めてみます.

「ここは,こうしたほうがいいよ.」などご指摘ありましたら,コメント欄などでご指導頂けると嬉しいです.

1. 1 変数の場合

簡単のために,1 変数の場合から考えていきましょう.$x$ 軸上の点 $\bar x$ 周りでの関数 $f(x)$ の値は,テイラー展開すると,

$$
f(\bar x + \Delta x) = f(\bar x) + \frac{df}{dx}(\bar x)\Delta x + \frac{1}{2}\frac{d^2f}{dx^2}(\bar x)\Delta x^2 + O(x^3)
$$

となります.($~O(\cdot)$ はランダウの O 記号です.) ここで,$3$ 次以上の項を無視すると,

$$
f(\bar x + \Delta x) = f(\bar x) + \frac{df}{dx}(\bar x)\Delta x + \frac{1}{2}\frac{d^2f}{dx^2}(\bar x)\Delta x^2
$$

となります.これを最大にすることを考えてみましょう.$\Delta x$ で微分して,$0$ とおくと以下の方程式を得ることができます.

$$
\frac{df}{dx}(\bar x) + \frac{d^2f}{dx^2}(\bar x)\Delta x = 0 \tag{1}
$$

ここで,$f^{\prime}(\bar x) = \cfrac{df}{dx}(\bar x),~~f^{\prime\prime}(\bar x) = \cfrac{d^2f}{dx^2}(\bar x)$ とおいて,$\Delta x$ について解くと,

$$
\Delta x = - \frac{f^{\prime}(\bar x)}{f^{\prime\prime}(\bar x)}
$$

となります.ここで,$\Delta x = x - \bar x$ であることから,

$$
x = \bar x - \frac{f^{\prime}(\bar x)}{f^{\prime\prime}(\bar x)}
$$

が得られます.$f(x)$ の 2 次近似 (2 次関数による近似) をして,その 2 次関数の極値を与えるような $x$ の値を求めているというのが,ここで行っている計算の意味です.次に,この極値を与えるような $x$ の値まわりで 2 次近似をして,さらに極値となる値を探します.これを繰り返すことで,徐々に $f(x)$ の極値に近づいていくというのがニュートン法のアイディアです.

まとめると,以下のような手続きになります.


ニュートン法 (1)
1. $x$ の初期値を与える
2. $\bar x \leftarrow x$ として,次のように更新する.
$$
x \leftarrow \bar x - \cfrac{f^{\prime}(\bar x)}{f^{\prime\prime}(\bar x)}
$$ 3. $~|~x - \bar x~| \leq \delta~$ かどうかを確認する.偽なら,step 2 に戻る.真なら x を返す.


1.1 例

参考書にあった解析的に解けて簡単な例でニュートン法を試してみましょう.以下のような関数 $f$ を考えます.

$$
f(x) = x^3 - 2 x^2 + x + 3

$$

これの 1 階導関数と 2 階導関数は,

$$
\frac{df}{dx}(x) = 3x^2 - 4 x + 1\\
\frac{d^2f}{dx^2}(x) = 6x - 4
$$

です.では,2 次近似を考えましょう.

$$
f(\bar x + \Delta x) = f(\bar x) + \frac{df}{dx}(\bar x)\Delta x + \frac{d^2f}{dx^2}(\bar x)\Delta x^2
$$

$\Delta x = x - \bar x$ なので,

$$
f(x) = f(\bar x) + \frac{df}{dx}(\bar x)(x - \bar x) + \frac{d^2f}{dx^2}(\bar x)(x - \bar x)^2
$$

となります.では,$f$ とその $\bar x = 2$ 周りの 2 次近似を同平面にプロットしてみます!

f = lambda x: x**3 - 2*x**2 + x + 3
d1f = lambda x: 3*x**2 - 4*x + 1
d2f = lambda x: 6*x - 4
f2 = lambda bar_x, x: f(bar_x) + d1f(bar_x)*(x - bar_x) + d2f(bar_x)*(x - bar_x)**2

x = np.linspace(-10, 10, num=100)
plt.plot(x, f(x), label="original")
plt.plot(x, f2(2, x), label="2nd order approximation")
plt.legend()
plt.xlim((-10, 10))

赤の曲線が元の関数 $f$ で,青の 2 次関数が近似した曲線になります.それでは,ニュートン法を Python で実装して停留点を与える解を求めます.本記事では,scipy.optimize など便利な API は敢えて使わずに作ります.

newton_1v.py
class Newton(object):
    def __init__(self, f, d1f, d2f, f2):
        self.f = f
        self.d1f = d1f
        self.d2f = d2f
        self.f2 = f2

    def _update(self, bar_x):
        d1f = self.d1f
        d2f = self.d2f
        return bar_x - d1f(bar_x)/d2f(bar_x)

    def solve(self, init_x, n_iter, tol):
        bar_x = init_x
        for i in range(n_iter):
            x = self._update(bar_x)
            error = abs(x - bar_x)
            print("|Δx| = {0:.2f}, x = {1:.2f}".format(error, x))
            bar_x = x
            if error < tol:
                break
        return x

def _main():
    f = lambda x: x**3 - 2*x**2 + x + 3
    d1f = lambda x: 3*x**2 - 4*x + 1
    d2f = lambda x: 6*x - 4
    f2 = lambda bar_x, x: f(bar_x) + d1f(bar_x)*(x - bar_x) + d2f(bar_x)*(x - bar_x)**2

    newton = Newton(f=f, d1f=d1f, d2f=d2f, f2=f2)
    res = newton.solve(init_x=10, n_iter=100, tol=0.01)
    print("Solution is {0:.2f}".format(res))

if __name__ == "__main__":
    _main()

これを実行すると,以下のような出力を得ることができます.

|Δx| = 4.64, x = 5.36
|Δx| = 2.30, x = 3.06
|Δx| = 1.10, x = 1.96
|Δx| = 0.47, x = 1.48
|Δx| = 0.14, x = 1.35
|Δx| = 0.01, x = 1.33
|Δx| = 0.00, x = 1.33
Solution is 1.33

プロットすると,以下のようになります.0 に向かって単調減少していますね.実際,得られた解は 1.33 で,これは解析解の $4/3$ に相当するものです.初期値を変えることで $0$ を得ることもできます.ただし,逆を返せば「初期値によって得られる解が異なる」ということでもあります.

error.png

ここまでのコードは,Gist に上げてありますので詳しくは こちら を御覧ください.

2. 多変数の場合

次に,多変数の場合について考えていきましょう.点 $[\bar x_1,~\cdots,~\bar x_n]$ の近くの点 $[\bar x_1 + \Delta x_1,~\cdots,~\bar x_n + \Delta x_n]$ での関数 $f(\bar x_1,~\cdots,~\bar x_n)$ の値はテイラー展開すると,
$$
f(\bar x_1 + \Delta x_1,~\cdots,~\bar x_n + \Delta x_n) = \bar f + \sum^{n}_{i=1} \frac{\partial \bar f}{\partial x_{i}}\Delta x_i + \frac{1}{2!} \sum^{n}_{i,~j=1} \frac{\partial^2 \bar f}{\partial x_i\partial x_j}\Delta x_i \Delta x_j + \cdots
$$

となります.この展開の導出は本記事の主旨から漏れてしまうので行いませんが,樋口三郎先生の講義 2 がとてもわかりやすいので参考にされると良いと思われます.

さて,テイラー展開に戻ります.先ほどのテイラー展開で用いた $\bar f$ は点 $[\bar x_1,~\cdots,~\bar x_n]$ における関数 $f$ の値になります.1 変数の場合と同様,3 次以降の項は無視して,$\Delta x_i$ で微分した後に $0$ とおきます.

$$
\frac{\partial \bar f}{\partial x_i} + \sum_{j=1}^{n}\frac{\partial^2 \bar f}{\partial x_i\partial x_j}\Delta x_j = 0
$$

これは他の変数に関しても同様であるので,ヘッセ行列 3 を考えると,

$$
\nabla \bar f + \boldsymbol{\bar H}\Delta\boldsymbol{x} = \boldsymbol{0}
$$

です.$\nabla$ はナブラという記号で,例えば 2 変数のスカラー関数を例に取ると,以下のように 1 階偏微分を要素に持つようなベクトルを作ります.これを勾配と言ったります.

$$
\nabla f(x_0,~x_1) = \begin{bmatrix}
\cfrac{\partial \bar f}{\partial x_0} \\
\cfrac{\partial \bar f}{\partial x_1}
\end{bmatrix}
$$

式 (1) と比較すると,勾配ベクトルが 1 階導関数,ヘッセ行列が 2 階導関数に対応しているのがわかるかと思います.これを $\Delta \boldsymbol{x} = \boldsymbol{x} - \boldsymbol{\bar x}$ に注意して,$\boldsymbol{x}$ に関して解きましょう.

$$
\boldsymbol{x} = \boldsymbol{\bar x} - \boldsymbol{\bar H}^{-1}\nabla \bar f
$$

以上から,このときのニュートン法は以下のような手続きとなります.


ニュートン法 (2)
1. $\boldsymbol{x}$ の初期値を与える
2. 勾配 $\nabla f$ とヘッセ行列 $\boldsymbol{H}$ の $\boldsymbol{\bar x}$ における値を計算.
3. 次の線形方程式の解を計算.
$$
\boldsymbol{H}\Delta\boldsymbol{x} = - \nabla f
$$
4. $\boldsymbol{x}$ を更新.
$$
\boldsymbol{x}\leftarrow \boldsymbol{x} + \Delta \boldsymbol{x}
$$
5. $||\Delta\boldsymbol{x}|| \leq \delta$ なら $\boldsymbol{x}$ を返す.偽なら Step2 に戻る.


2.1 例

では,参考書 1 でも取り上げられている次のような関数を考えましょう.

$$
f(\boldsymbol{x}) = x_{0}^{3} + x_{1}^{3} - 9x_{0}x_{1} + 27
$$

これの $\boldsymbol{\bar x} = [\bar x_0, \bar x_1]$ 周りでの 2 次近似は,

$$
f(\boldsymbol{x}) = \bar f + \sum^{1}_{i=0}\frac{\partial \bar f}{\partial x_i}(x_i - \bar x_i) + \frac{1}{2!}\sum^{1}_{i,~j=0}\frac{\partial^2 \bar f}{\partial x_i\partial x_j}(x_i - \bar x_i)(x_j - \bar x_j)
$$

で,勾配 $\nabla f$ とヘッセ行列 $\boldsymbol{H}$ は,

$$
\nabla f = \begin{bmatrix}
\partial \bar f/\partial x_0 \\
\partial \bar f/\partial x_1
\end{bmatrix} = \begin{bmatrix}
3 \bar x_0^2 - 9 \bar x_1 \\
3 \bar x_1^2 - 9 \bar x_0
\end{bmatrix}\\
\boldsymbol{H}= \begin{bmatrix}
\partial^2 \bar f/\partial x_0^2 & \partial^2 \bar f/\partial x_0\partial x_1 \\
\partial^2 \bar f/\partial x_1\partial x_0 & \partial^2 \bar f/\partial x_1^2
\end{bmatrix} = \begin{bmatrix}
6\bar x_0 & -9 \\
-9 & 6\bar x_1
\end{bmatrix}
$$

となります.では,Python で実装してみましょう.

newton_mult.py
import numpy as np
import matplotlib.pyplot as plt

plt.style.use("ggplot")
plt.rcParams["font.size"] = 13
plt.rcParams["figure.figsize"] = 16, 12

class MultiNewton(object):
    def __init__(self, f, dx0_f, dx1_f, grad_f, hesse):
        self.f = f
        self.dx0_f = dx0_f
        self.dx1_f = dx1_f
        self.grad_f = grad_f
        self.hesse = hesse

    def _compute_dx(self, bar_x):
        grad = self.grad_f(bar_x)
        hesse = self.hesse(bar_x)
        dx = np.linalg.solve(hesse, -grad)
        return dx

    def solve(self, init_x, n_iter=100, tol=0.01, step_width=1.0):
        self.hist = np.zeros(n_iter)
        bar_x = init_x
        for i in range(n_iter):
            dx = self._compute_dx(bar_x)
            # update
            x = bar_x + dx*step_width
            print("x = [{0:.2f} {1:.2f}]".format(x[0], x[1]))

            bar_x = x
            norm_dx = np.linalg.norm(dx)
            self.hist[i] += norm_dx
            if norm_dx < tol:
                self.hist = self.hist[:i]
                break
        return x

def _main():
    f = lambda x: x[0]**3 + x[1]**3 - 9*x[0]*x[1] + 27
    dx0_f = lambda x: 3*x[0]**2 - 9*x[1]
    dx1_f = lambda x: 3*x[1]**2 - 9*x[0]
    grad_f = lambda x: np.array([dx0_f(x), dx1_f(x)])
    hesse = lambda x: np.array([[6*x[0], -9],[-9, 6*x[1]]])

    init_x = np.array([10, 8])
    solver = MultiNewton(f, dx0_f, dx1_f, grad_f, hesse)
    res = solver.solve(init_x=init_x, n_iter=100, step_width=1.0)
    print("Solution is x = [{0:.2f} {1:.2f}]".format(res[0], res[1]))

    errors = solver.hist
    epochs = np.arange(0, errors.shape[0])

    plt.plot(epochs, errors)
    plt.tight_layout()
    plt.savefig('error_mult.png')

if __name__ == "__main__":
    _main()

これを実行すると,以下のような出力を得ることができます.解は $[3, 3]$ となり解析解の 1 つを得ることができています.

x = [5.76 5.08]
x = [3.84 3.67]
x = [3.14 3.12]
x = [3.01 3.00]
x = [3.00 3.00]
Solution is x = [3.00 3.00]

最後に,どのように収束していくのか $||\Delta x||$ の推移のグラフを示して終わります.

error_mult.png

おわりに

以上,1 変数にはじまって例が 2 変数ではありますが多変数の場合のニュートン法について概説しました.機械学習のアルゴリズムを勉強していると,何かしらの基準値を定義してその関数の最小値 (または最大値) を与えるパラメータを求める場面に多く当たると思われます.このときに,こういったニュートン法などの最適化手法を用いるわけですね.数式をバーっと眺めると難しい印象がありますが,一から読んでみると大学数学あたりの知識があれば意外と追えることもわかります.

近年では,scikit-learn などを使えばすぐに機械学習のアルゴリズムを試せますが,中身がどうなっているのかを理解することも大事ですよね.と言ったところで,この記事を終わりにします.

(気が向いたら 準ニュートン法なども書ければと.)

References