共役勾配法によるコスト関数最適化 - 人工知能に関する断創録

今回もCourseraの機械学習ネタ。Courseraの講義ではロジスティック回帰やニューラルネットのパラメータ学習に共役勾配法（conjugate gradient method: CG法）やBFGSアルゴリズム（Broyden–Fletcher–Goldfarb–Shanno algorithm）*1 を使っていました。よく使われる勾配降下法（gradient descent method）より強力な最適化アルゴリズムとのこと。勾配降下法に比べてよい点として

学習率を手動で設定する必要がない
高速

という2点が挙げられていました。自分で試してみるとよくわかるけど、学習率の決定は地味に面倒くさいので非常に助かります。デメリットとしては実装が複雑なことが挙げられていました。でも普通は数値計算のプロが書いたライブラリ使うからほとんどメリットしかないね！CourseraではMatlabまたはOctaveのfminunc()を使っていましたが、今回はPythonのscipy.optimizeの関数を使って書き直してみます。

コスト関数と導関数

関数最適化に関するライブラリはscipy.optimizeモジュールにまとまっていました。いろいろありますが、今回は共役勾配法のfmin_cg()とBFGSアルゴリズムのfmin_bfgs()を使ってみます。この関数に与えるのは最適化したいコスト関数とその導関数です。コスト関数のパラメータがn個ある場合は、それぞれのパラメータで偏微分するのでn個の導関数のリストになります。

$\displaystyle J(\theta)$
$\displaystyle \frac{\partial}{\partial \theta_j} J(\theta) \hspace{25pt} (j=1, \cdots, n)$

CG法やBFGSアルゴリズムでは、上のn+1個の関数とパラメータの初期値を与えると $J(\theta)$ を最小化するパラメータ $\theta$ を求めてくれます。つまり、 $min_\theta J(\theta)$ です。いきなりロジスティック回帰やニューラルネットのコスト関数の最適化は敷居が高いので答えがわかっている簡単なコスト関数を最適化してパラメータを求めてみます。では早速やってみます。

（例1）

パラメータ $\theta$ とコスト関数 $J(\theta)$ を下のように定義します。

$\displaystyle \theta = [ \theta_1 \theta_2 ]^T$
$\displaystyle J(\theta) = (\theta_1 - 5)^2 + (\theta_2 - 5)^2$

それぞれのパラメータでコスト関数を偏微分すると下の二つの導関数が得られます。

$\displaystyle \frac{\partial}{\partial \theta_1} J(\theta) = 2 (\theta_1 - 5)$
$\displaystyle \frac{\partial}{\partial \theta_2} J(\theta) = 2 (\theta_2 - 5)$

それぞれの導関数を＝0としたときのパラメータが最適解なので、答えは $\theta_1 = 5$ 、 $\theta_2 = 5$ とすぐにわかります。では、これをfmin_cg()とfmin_bfgs()で解けるか確かめてみます。

#coding: utf-8
import numpy as np
from scipy import optimize

def J(theta, *args):
    """最小化を目指すコスト関数を返す"""
    theta1, theta2 = theta
    return (theta1 - 5) ** 2 + (theta2 - 5) ** 2

def gradient(theta, *args):
    """コスト関数の偏微分を返す
    各パラメータで偏微分した関数リストを返す"""
    theta1, theta2 = theta
    # Jをtheta1で偏微分した関数
    gt1 = 2 * (theta1 - 5)
    # Jをtheta2で偏微分した関数
    gt2 = 2 * (theta2 - 5)
    return np.array((gt1, gt2))

if __name__ == "__main__":
    # パラメータの初期値
    initial_theta = np.array([0, 0])

    # Conjugate Gradientによるコスト関数最小化
    theta = optimize.fmin_cg(J, initial_theta, fprime=gradient)
    print "conjugate gradient:", theta

    # BFGSによるコスト関数最小化
    theta = optimize.fmin_bfgs(J, initial_theta, fprime=gradient)
    print "BFGS:", theta

実行結果は、

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 2
         Function evaluations: 5
         Gradient evaluations: 5
conjugate gradient: [ 4.99999995  4.99999995]
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 2
         Function evaluations: 4
         Gradient evaluations: 4
BFGS: [ 5.  5.]

どちらも、 $\theta_1 = 5$ 、 $\theta_2 = 5$ を返していることがわかります。

（例2）

次はもう少し複雑な例を。scipy.optimizeの例に載っていたサンプルです。J()とgradient()は最適化したいパラメータの関数にする必要があるみたいなので、それ以外のデータはまとめてargsに渡す必要があります。ロジスティック回帰だったらパラメータ推定に使用する教師データ (X, y) などです。今回は例題として関数の係数を引数としてargsに渡してみます。パラメータを (u, v) とし、コスト関数 $J(u, v)$ を下のように定義します。

$\displaystyle J(u, v) = a u^2 + b u v + c v^2 + d u + e v + f$

それぞれのパラメータでコスト関数を偏微分すると下の二つの導関数が得られます。

$\displaystyle \frac{\partial}{\partial u} J(u, v) = 2 a u + b v + d$
$\displaystyle \frac{\partial}{\partial v} J(u, v) = b u + 2 c v + e$

例として係数を a= 2, b = 3, c = 7, d = 8, e = 9, f = 10 とし、それぞれの導関数を＝0とした連立方程式を解くと答えは u = -1.8085106、v = -0.2553191 とすぐにわかります。では、これをfmin_cg()とfmin_bfgs()で解けるか確かめてみます。

#coding: utf-8
import numpy as np
from scipy import optimize

def J(param, *args):
    """最小化を目指すコスト関数を返す"""
    u, v = param
    # パラメータ以外のデータはargsを通して渡す
    a, b, c, d, e, f = args
    return a*u**2 + b*u*v + c*v**2 + d*u + e*v + f

def gradient(param, *args):
    """コスト関数の偏微分を返す
    各パラメータで偏微分した関数リストを返す"""
    u, v = param
    a, b, c, d, e, f = args
    gu = 2*a*u + b*v + d
    gv = b*u + 2*c*v + e
    return np.array((gu, gv))

if __name__ == "__main__":
    # パラメータの初期値
    initial_param = np.array([0, 0])
    args = (2, 3, 7, 8, 9, 10)

    # Conjugate Gradientによるコスト関数最小化
    param = optimize.fmin_cg(J, initial_param, fprime=gradient, args=args)
    print "conjugate gradient:", param

    # BFGSによるコスト関数最小化
    param = optimize.fmin_bfgs(J, initial_param, fprime=gradient, args=args)
    print "BFGS:", param

実行結果は、

Optimization terminated successfully.
         Current function value: 1.617021
         Iterations: 2
         Function evaluations: 5
         Gradient evaluations: 5
conjugate gradient: [-1.80851064 -0.25531915]
Optimization terminated successfully.
         Current function value: 1.617021
         Iterations: 2
         Function evaluations: 5
         Gradient evaluations: 5
BFGS: [-1.80851064 -0.25531915]

となり、手計算と同じになります。これはすごい・・・

最適化アルゴリズムの詳細を調べてもいいけれど、けっこう難しいし、あまり興味がわかない。ブラックボックスとして使いこなせれば十分かな。導関数を自分で計算しなくちゃいけないのがちょっと面倒くさい。Sympyとかで自動的に計算できないのかな？

あとでより複雑なコスト関数を持つロジスティック回帰やニューラルネットのパラメータ最適化にも応用してみたいと思います！

*1:Shannonではなく、Shannoさんなのか。勘違いしていた。