関数の最適化 ～勾配法～

バンドル調整に向けて，数学の復習をしておきます．参考文献はいつものごとく金谷先生の本です．

最適化とは

ある制約条件のもとで，ある関数の値を最大または最小にするような変数の値を求めること．

最大値・最小値を求めるためのトリック

ニュートン法，勾配法，レーベンバーグ・マーカート法等，非線形最適化をするアルゴリズムはたくさんありますが，基本となっている考え方は，，，

関数 f(x) の極値では，一次導関数 f'(x) = 0 となる．

です．結局アルゴリズムが違っても，f'(x) = 0 となるような x を探しているという点では同じです．なので，陰関数 f(x) = 0 を解くことができるツールであれば，最適化にも使えると思ってよいと思います．最適化の場合，解くべき陰関数が関数 f(x) ではなく，一次導関数 f'(x) となります．

一変数の場合の勾配法

まず導入として，一変数関数の場合の勾配法を考えます．基本となる考え方はシンプルで，上に凸の関数の場合，勾配がプラス（マイナス）の間はずっと x を増加（減少）させます．上に凸の関数であれば，この作業を続けていればどこかで必ず勾配がマイナス（プラス）になります．この間のどこかに勾配が０である場所があると考えられるので，これを求める最適解として採用します．

一変数の場合の勾配法アルゴリズム

最大値を求める場合の勾配法アルゴリズムは下記のとおりです．傾きの符号が変化するまで同一方向に進み，符号変化を検知した時点で進む距離（ステップ幅）を減らして反対方向に進む．という作業を繰り返し，最後には f'(x) = 0 を満たす変数を求めます．

多変数の場合の勾配法

少し複雑になりますが，多変数の場合も考え方は同じです．多変数の場合も曲面に対して接平面を考え，この接平面の傾きが一番大きい方向を基準に進んでいきます．で，この方向に進んでいく中で一番大きい値を探索します．最大値を見つけたら，その点に関して接平面を考え再び傾きが一番大きい方向を基準に進んでいき．．．．と繰り返せば最大値が求まります．

多変数の場合の勾配法アルゴリズム

一変数の勾配法の実験 with Python

ということで，勾配法を実装してみました．まずは一変数の場合の勾配法です．最大値を求めたい関数として，下記の関数を考えます．

f(x) = sin(x) + cos(x/2 -1)

一変数の勾配法による実験結果

下記実験結果ですが，緑色の点が反復の途中仮定，赤色の点が最終結果になります．最大値の回りで反復が終わっていることがわかります．

多変数の勾配法の実験 with Python

ついで，多変数の勾配法です．変数の数が多いと可視化が難しいので，二変数の関数に対して最大値を求めます．

f(x, y) = -(15*(x-6)**2 + 2*(x-6)*(y-4) + 3*(y-4)**2 + 4*(x-6) + 3*(y-4) - 15)

多変数の勾配法による実験結果

下記実験結果ですが，楕円のてっぺんが最大値になります．各点で等高線に垂直な方向に進み，最終的には最大値にたどり着いている様子がわかります．緑色の点が途中経過で，黄色の点が最終解です．

Pythonコード

下記，まだちょっとおかしいです．．．．収束しないときが多々あり．．．また時間ができたら見直します．

gradient_method.py

# -*- coding: utf-8 -*-

import numpy as np

def derivativeNd(func, xVec, d):
  dVec = np.zeros(xVec.shape)
  
  for i in range(xVec.shape[0]):
    xVecP = np.copy(xVec)
    xVecN = np.copy(xVec)
    
    xVecP[i] = xVec[i] + d/2
    xVecN[i] = xVec[i] - d/2

    dVec[i] = (func(xVecP) - func(xVecN)) / d

  return dVec

def differentiate(func, xVec, dirVec, d):

  dVec = np.zeros(xVec.shape)  
  dx = dirVec / np.linalg.norm(dirVec) * d
  
  for i in range(xVec.shape[0]):
    xVecP = np.copy(xVec)
    xVecN = np.copy(xVec)
    
    xVecP[i] = xVec[i] + dx[i] / 2
    xVecN[i] = xVec[i] - dx[i] / 2
  
    dVec[i] = (func(xVecP) - func(xVecN)) / d  

  return dVec
  
def hill_climbNd(func, xVecInit, stepInit, ipu):
  
  xVec = np.copy(xVecInit)
  xVecOld = np.copy(xVecInit)
  step = stepInit
  
  history = np.zeros((0, xVec.shape[0]))
  history = np.append(history, np.array([xVec]), axis=0)
  
  while True:
    xVecOld = np.copy(xVec)
    dirVec = derivativeNd(func, xVec, 0.000001)
    xVec, tmp = hill_climb_to_dir(func, xVec, dirVec, step, ipu)
    history = np.append(history, np.array([xVec]), axis=0)
    
    if np.linalg.norm(xVec - xVecOld) <= ipu:
      break
  
  return xVec, history

def hill_climb_to_dir(func, xVecInit, dirVec, stepInit, ipu):

  x = xVecInit
  step = stepInit
  d = 0.00001
  normDirVec = dirVec / np.linalg.norm(dirVec)
  history = [x]
  MINSTEP = 1.0e-15  

  while True:
    # Step 2 : Setup Step Length
    dfVec = differentiate(func, x, normDirVec, min(d, abs(step)/1000))
    
    if np.linalg.norm(dfVec) <= ipu:
      break
    
    sign = 1 if np.dot(dfVec, normDirVec) > 0 else -1   
    xb = x
    xa = x + sign * abs(step) * normDirVec
    
    # Step 3 : Find Another Side of Mountain.
    if func(xb) < func(xa):
      while True:
        step = 2*step
        xb = xa
        xa = xb + sign * abs(step) * normDirVec
        
        if func(xb) >= func(xa):
          break
      x = xb
      step = max(step/2, MINSTEP)
        
    # Step 4 : Overrun and Reset.
    else:
      while True:
        step = max(step/2, MINSTEP)
        xa = xa - sign * abs(step) * normDirVec
        
        if func(xb) <= func(xa):
          break
          
      x = xa
      step = 2*step
    
    # Step 5 : Convergence Check
    dfVec = differentiate(func, x, normDirVec, min(d, abs(step)/1000))   
    history.append(x)
    
    if np.linalg.norm(dfVec) < ipu:
      break
  
  return x, np.array(history)

test_gradient_method.py

# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt 
import numpy as np
import math

from optimization.gradient_method import hill_climb_to_dir
from optimization.gradient_method import hill_climbNd

def test_func1d(xVec):
  return math.sin(xVec[0]) + np.cos(xVec[0]/2 - 1)

def test_func2d(xVec):
  x = xVec[0]
  y = xVec[1]
  return (-(15*(x-6)**2 + 2*(x-6)*(y-4) + 3*(y-4)**2 + 4*(x-6) + 3*(y-4) - 15))
  
def test_hill_climb1d():
  
  x = np.linspace(-5, 5, 101)
  y = np.zeros(x.shape[0])
  for i in range(x.shape[0]):
    y[i] = test_func1d(np.array([x[i]]))
    
  ans, hist = hill_climb_to_dir(test_func1d, np.array([-2.0]),\
                                np.array([1.0]), 0.1, 0.00001)
  
  xVecInit = np.array([-2.0], dtype='float64')
  stepWidth = 1.0
  termCondition = 0.01
  xVec, history = \
    hill_climbNd(test_func1d, xVecInit, stepWidth, termCondition)
  
  fig = plt.figure()
  ax = fig.add_subplot(1, 1, 1)
  ax.plot(x, y)  
  
  for i in range(hist.shape[0]):
    x = hist[i]
    y = test_func1d(x)
    if i == hist.shape[0] - 1:
      size = 50
      color = 'red'
    else:
      size = 10
      color = 'green'
      
    ax.scatter(x, y, s=size, c=color, edgecolors=color)
    
  plt.show()

  
def test_hill_climb2d():
  
  x1 = np.linspace(-10, 10, 101)
  x2 = np.linspace(-10, 10, 101)
  X1, X2 = np.meshgrid(x1, x2)
  X = np.c_[np.ravel(X1), np.ravel(X2)]
  
  Z = np.empty((0, 1))
  for i in range(X.shape[0]):
    z = test_func2d(X[i])
    Z = np.append(Z, [[z]], axis=0)
  Z = np.reshape(Z, (X1.shape[0], X2.shape[0]))

  xVecInit = np.array([-3.0, -3.0], dtype='float64')
  stepWidth = 1.0
  termCondition = 0.01
  xVec, history = \
    hill_climbNd(test_func2d, xVecInit, stepWidth, termCondition)
  
  fig = plt.figure()
  ax = fig.add_subplot(1, 1, 1)

  values = set()
  for i in range(history.shape[0]):
    xVec = history[i]
    values.add(test_func2d(xVec))
  
  cont = ax.contourf(X1, X2, Z)
  for i in range(cont.levels.shape[0]):
    values.add(cont.levels[i])
  
  levels = list(values)
  levels.sort()
  ax.contourf(X1, X2, Z, levels)

  for i in range(history.shape[0]):
    xVec = history[i]

    if i == history.shape[0] - 1:
      size = 50
      color = 'yellow'
    else:
      size = 10
      color = 'green'      
    ax.scatter(xVec[0], xVec[1], s=size, c=color, edgecolors=color)
    
  plt.show()
  
if __name__ == "__main__":
  print("Gradient Method Test")
  
  #test_hill_climb1d()  
  test_hill_climb2d()

参考文献

これなら分かる最適化数学―基礎原理から計算手法まで