-Python- 勾配降下法

下式のような制約条件のない最適化問題を考えます.

Minimize  5x^2 - 6xy + 3y^2 + 6x -6y

ここで,[TeX: f(x, y) = 5x^2 - 6xy + 3y^2 + 6x -6y] とし,[TeXf(x, y) = k] を満たす点の集合を考えます.
の勾配は

f:id:HidehikoMURAO:20181106175937p:plain
となり,点 (x_0, y_0)においては,その点を通る等高線の接戦に垂直方向で,kが大きくなる方向を向いたベクトルが ∇ f (x_0, y_0)になります.kを小さくするには,-∇ fの方向に進めば良いことになります.
 
あるパラメータαと初期点x_0 = (x_0, y_0)を最初に決めておいて,現在の地点x_kからx_k+1 = x_k - α∇f(x_k) に移動するということを繰り返して行けば,目的関数の値は小さくなるはずです.局地においては∇ f (x, y) = 0 となるので,そのような点が見つかれば計算が終了します.
 
上記の2変数関数の場合を拡張して,一般にn変数関数で以下のような同様のアルゴリズム(勾配降下法または最急降下法)を考えることが可能です.
  • パラメータ αεが入力値として与えられる.
  • 初期点 x_k を決める.
  • k を0から1づつ増やしながら以下を繰り返す.
    • || ∇ f (x_k) || ≦ ε であれば終了.
    • x_k+1 = x_k - α∇ f (x_k)を計算する.

このアルゴリズムにおいて,パラメータαは最適値の検索にどのくらい大きく移動するかを表します.αとしてどのような値が良いかは,問題によって実験的に求められます.この値はあまり大きすぎると計算が発散してしまいますし,小さすぎるとなかなか収束しません.
パラメータεは計算の終了条件を決めていて,小さくすればするほど正確な最適解が得られますが,計算に時間がかかるようになります(例えば,判定条件を零に完全一致することを終了条件とすると計算が不安定になり,無限ループに陥ってしまったりする).

以下にアルゴリズムの例を示します.

 
# Import Module
import numpy as np
 
# Definition of Class
class GradientDescent: 
    def __init__(self, f, df, alpha = 0.01, eps = 1e-6): 
        self.f = f
        self.df = df
        self.alpha = alpha
        self.eps = eps
        self.path = None
    
    def solve(self, init): 
        x = init
        path = []
        grad = self.df(x)
        path.append(x)
        while (grad**2).sum() > self.eps**2
            x = x - self.alpha * grad
            grad = self.df(x)
            path.append(x)
        self.path_ = np.array(path)
        self.x_ = x

 

        self.opt_ = self.f(x)
  

このクラスは変数の数が一般の場合に解けるようになっています.コンストラクタの引数の f とdf には,最小化したい関数と導関数を指定します.n変数関数の最適化をする場合,fはn次元ベクトル(長さがnの1次元配列)を引数にとり,浮動小数点型の戻り値を返します.dfはn次元ベクトルを引数にとり,戻り値もn次元ベクトルを返します.そのようなfとdfを引数として与えることで,一般の場合に使用することができます.

デフォルト値付きの引数 alpha は,アルゴリズムの中のαで,探索時の移動の大きさを表します.引数epsはアルゴリズムの終了条件の基準を表します. \nabla f のL2ノルムがeps以下になった時点で終了します.メソッドsolveでは引数として初期値(アルゴリズム中の x_0)をとり,計算結果としての最適解はデータ属性x_に保存され,最適値はデータ属性opt_に保存されます.

 

実際に  f(x, y) = 5x^2 - 6xy + 3y^2 + 6x -6y の最小値を計算します.まずこれを偏微分すると 

 f:id:HidehikoMURAO:20181106175552p:plain

となります.これをもとに以下のように実装します.最適値の計算のみではなく,計算途中の様子も可視化できるようになっています.

# Import Module
import numpy as np
import matplotlib.pyplot as plt
import gd
 
# Definition of Function
def f(xx): 
    x = xx[0]
    y = xx[1]
    return 5 * x**2 - 6*x*y + 3*y**2 + 6*x - 6*y
 
def df(xx): 
    x = xx[0]
    y = xx[1]
    return np.array([10*x - 6*y + 6, -6 * x + 6*y - 6])
 
# Solving problem
algo = gd.GradientDescent(f, df) 
initial = np.array([1, 2])
algo.solve(initial)
 
# Display results
print(algo.x_)
print(algo.opt_)
 
# Drawing Graph
plt.scatter(initial[0], initial[1], color = "k", marker = "o")
plt.plot(algo.path_[:, 0], algo.path_[:, 1], color = "k"
         linewidth=1.5)
xs = np.linspace(-2, 2, 300)
ys = np.linspace(-2, 2, 300)
xmesh, ymesh = np.meshgrid(xs, ys)
xx = np.r_[xmesh.reshape(1, -1), ymesh.reshape(1, -1)]
levels = [-3, -2.9, -2.8, -2.6, -2.4
          -2.2, -2, -1, 0, 1, 2, 3, 4]
plt.contour(xs, ys, f(xx).reshape(xmesh.shape), 
            levels = levels, colors = "k", linestyles = "dotted")
 
plt.show()