JavaScriptを有効にしてください

直線探索を使った最急降下法をPythonで実装

 ·   9 min read

※記事内に商品プロモーションを含むことがあります。

はじめに

本記事では、最急降下法と、Armijo条件と呼ばれる直線探索手法について簡単に解説する。

数理工学社の「工学基礎 最適化とその応用」を読んだので、4章「非線形計画法I(無制約最小化問題)」から、直線探索を使った最急降下法をPythonで実装した。

最急降下法とは勾配を用いる最適化手法の1つであり、最急降下法で解の更新幅を求める手法を直線探索という。
本記事では、これらについて簡単に解説し、Pythonで実装したコードを示す。
最後に、実装したコードを使って、簡単な2変数関数を最小化する。

また、以降ではライブラリを次のようにインポートしていることを前提とする。

1
2
import numpy as np
import matplotlib.pyplot as plt

バージョンは以下の通り。

ソフトウェア バージョン
python 3.7.4
numpy 1.16.4
matplotlib 3.1.0

最急降下法

最急降下法 (gradient descent, またはsteepest descent) は、数値最適化手法の1つであり、関数の勾配の方向に解を更新して、最適解に収束させようとするものである。

まず、nn次元ベクトルxRn\boldsymbol{x} \in \mathbb{R}^nの関数f(x)Rf(\boldsymbol{x}) \in \mathbb{R}を最小化する問題を考える。
f(x)f(\boldsymbol{x})は微分可能であり、勾配ベクトルは次式で表されるとする。
f(x)=[fx1,,fxn]TRn \nabla f(\boldsymbol{x}) = \left[ \frac{\partial f}{\partial x_1}, …, \frac{\partial f}{\partial x_n} \right]^\mathrm{T} \in \mathbb{R}^n

次に、最急降下法のアルゴリズムについて述べる。

  1. 反復回数をk=0 k=0として、解の初期値xk\boldsymbol{x_k}を与える。
  2. 終了条件を満たせば最適解をxk\boldsymbol{x_k}として終了する。そうでなければ、xk\boldsymbol{x_k}における勾配f(xk) \nabla f(\boldsymbol{x_k})を求める。
  3. 解の探索方向をdk=f(xk) \boldsymbol{d_k} = -\nabla f(\boldsymbol{x_k})とする。
  4. 探索方向のステップ幅αk \alpha_kを直線探索により求め、新たな解をxk+1=xk+αkdk \boldsymbol{x_{k+1}} = \boldsymbol{x_k} + \alpha_k \boldsymbol{d_k}と更新する。
  5. kk+1 k \leftarrow k+1として、2に戻る。

下図は2変数関数における勾配降下法のイメージである。
等高線は関数の値を示し、(x1,x2)=(0,0) (x_1, x_2)=(0, 0)で最小値をとる。
解候補の点xk\boldsymbol{x_k}における探索方向dk \boldsymbol{d_k}は、勾配が最も急な方向(等高線と垂直な方向)となる。

gradient-descent

最急降下法には、ステップ幅αk \alpha_kが大きすぎると解が発散し、小さすぎると収束が遅くなる問題がある。
この問題に対処できるようにαk \alpha_kを決める手法が、次節の直線探索である。

直線探索

直線探索 (line search) は、最急降下法や準ニュートン法などの最適化手法において得られた解の探索方向に対して、適切な更新幅を求める手法である。
解を大域的に収束させることが求められる。
実用的な直線探索手法としてArmijo(アルミホ)条件とWolf(ウルフ)条件があり、本記事では前者のみ扱う。

Armijo条件とは、0<ξ<1 0<\xi<1である定数ξ\xiに対して、
f(xk+αdk)f(xk)+ξαf(xk)Tdk f(\boldsymbol{x_k}+\alpha \boldsymbol{d_k}) \leq f(\boldsymbol{x_k}) + \xi \alpha \nabla f(\boldsymbol{x_k})^\mathrm{T} \boldsymbol{d_k}
を満たすα>0 \alpha > 0を選ぶものである。

下図にArmijo条件を幾何的に示す。
関数y=f(xk+αdk) y=f(\boldsymbol{x_k}+\alpha \boldsymbol{d_k})dk \boldsymbol{d_k}方向の目的関数である。
Armijo条件を満たすα\alphaとは、目的関数y=f(xk+αdk) y=f(\boldsymbol{x_k}+\alpha \boldsymbol{d_k})が、直線f(xk)+ξαf(xk)Tdk f(\boldsymbol{x_k}) + \xi \alpha \nabla f(\boldsymbol{x_k})^\mathrm{T} \boldsymbol{d_k}よりも小さい区間にある。
なお、ξ=0 \xi=0のとき、直線の値はf(xk) f(\boldsymbol{x_k})と等しくなるので解の収束が保障されなくなる。
また、ξ=1 \xi=1のとき、直線は接線y=f(xk)+αf(xk)Tdk y=f(\boldsymbol{x_k}) + \alpha \nabla f(\boldsymbol{x_k})^\mathrm{T} \boldsymbol{d_k}と一致するのでα=0 \alpha=0となり、解の更新幅が0になる。

line-search-armijo

Armijo条件を満たすステップ幅αk \alpha_kを得るアルゴリズムは、以下の通り。

  1. パラメータ0<ξ<1,0<τ<1 0<\xi<1, 0<\tau<1を与える。
  2. 更新幅の初期値をα=1 \alpha=1とおく。
  3. α \alphaがArmijo条件を満たすならば終了。
  4. 満たさなければατα \alpha \leftarrow \tau \alphaと更新して3に戻る。

上記のアルゴリズムでは、初めにある程度大きなステップ幅を与え、Armijo条件を満たすまで徐々に小さくしている。そのため、有限回の反復でArmijo条件を満たすステップ幅を得られることが保証される。

Pythonによる実装

Armijo条件を用いた最急降下法をPythonで実装した。
funに目的関数、derに勾配を返す関数を与える。また、minimizeメソッドにxを初期値として与えると、最急降下法を実行する。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
class GradientDescent:
    def __init__(self, fun, der, xi=0.3, tau=0.9, tol=1e-6, ite_max=2000):
        self.fun = fun         # 目的関数
        self.der = der         # 関数の勾配
        self.xi  = xi          # Armijo条件の定数
        self.tau = tau         # 方向微係数の学習率
        self.tol = tol         # 勾配ベクトルのL2ノルムがこの値より小さくなると計算を停止
        self.path = None       # 解の点列
        self.ite_max = ite_max # 最大反復回数
        
    def minimize(self, x):
        path = [x]
        
        for i in range(self.ite_max):
            grad = self.der(x)
            
            if np.linalg.norm(grad, ord=2)<self.tol:
                break
            else:
                beta = 1
                
                while self.fun(x - beta*grad) > (self.fun(x) - self.xi*beta*np.dot(grad, grad)):
                    # Armijo条件を満たすまでループする
                    beta = self.tau*beta
                    
                x = x - beta * grad
                path.append(x)
        
        self.opt_x = x                # 最適解
        self.opt_result = self.fun(x) # 関数の最小値
        self.path = np.array(path)    # 探索解の推移

最適化の例

次式の2変数関数を最小化する例を示す。

f(x)=2x12+x22+x1x2,(x=[x1,x2]T) f(\boldsymbol{x}) = 2x_1^2 + x_2^2 + x_1 x_2, (\boldsymbol{x}=[ x_1, x_2 ]^\mathrm{T})

この関数は、(x1,x2)=(0,0) (x_1, x_2)=(0, 0)で最小値0をとる。
また、勾配ベクトルは次式で与えられる。

f(x)=[fx1,fx2]T=[4x1+x2,x1+2x2]T \nabla f(\boldsymbol{x}) = \left[ \frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2} \right]^\mathrm{T} = [4x_1 + x_2, x_1 + 2x_2 ]^\mathrm{T}

関数値と勾配を求める関数を、それぞれf(x), f_der(x)として定義する。

1
2
3
4
5
def f(x):
    return 2*x[0]**2 + x[1]**2 + x[0]*x[1]

def f_der(x):
    return np.array([4*x[0] + x[1], x[0] + 2*x[1]])

関数の等高線をプロットする。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
x1 = np.linspace(-2, 2, 21)
x2 = np.linspace(-2, 2, 21)
x1_mesh, x2_mesh = np.meshgrid(x1, x2)
z = f(np.array((x1_mesh, x2_mesh)))

fig, ax = plt.subplots(figsize=(6, 6))
ax.contour(x1, x2, z, levels=np.logspace(-0.3, 1.2, 10))
ax.set_xlabel("x1")
ax.set_ylabel("x2")
ax.set_aspect('equal')
plt.show()

line-search-function

また、勾配をプロットする。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
x1 = np.linspace(-2, 2, 21)
x2 = np.linspace(-2, 2, 21)
x1_mesh, x2_mesh = np.meshgrid(x1, x2)

grad = f_der(np.array((x1_mesh, x2_mesh)))
U = grad[0] # x1方向の勾配
V = grad[1] # x2方向の勾配

fig, ax = plt.subplots(figsize=(6, 6))
ax.quiver(x1, x2, U, V, color='blue')
ax.set_aspect('equal')
plt.show()

line-search-gradient

前節のGradientDescentクラスを用いて、最急降下法を実行する。
初期値は(x1,x2)=(1.5,1.5) (x_1, x_2)=(1.5, 1.5)とする。

1
2
3
gd = GradientDescent(f, f_der)
init = np.array([1.5, 1.5])
gd.minimize(init)

解の推移をプロットする。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
path = gd.path

x1 = np.linspace(-2, 2, 21)
x2 = np.linspace(-2, 2, 21)
x1_mesh, x2_mesh = np.meshgrid(x1, x2)
z = f(np.array((x1_mesh, x2_mesh)))

fig, ax = plt.subplots(figsize=(6, 6))
ax.contour(x1, x2, z, levels=np.logspace(-0.3, 1.2, 10))
ax.plot(path[:,0], path[:,1], marker="o")
ax.set_xlabel("x1")
ax.set_ylabel("x2")
ax.set_aspect('equal')
plt.show()

line-search-path

発散せずに最適解へ収束しており、Armijo条件を用いた直線探索がうまく機能していることが分かる。

参考

以下のページにステップ幅を固定した最急降下法のコードがあり、実装のベースにさせて頂いた。
最急降下法の概要 - Qiita

直線探索の詳細なアルゴリスムについては、以下の書籍を参考とした。

シェアする

Helve
WRITTEN BY
Helve
関西在住、電機メーカ勤務のエンジニア。X(旧Twitter)で新着記事を配信中です

サイト内検索