Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2017-05-02

MikuHatsune2017-05-02

勾配法による最適化

ある関数があって、最小化もしくは最大化したい。これらは正負を入れ替えればよいでの最小化を考える。

 

いま、適当にf(x,y)=x(x-1)(x-2)+y(y-2)(y-3)(y+7)-8 という関数があったとする。これの最小値とそのときの(x,y)の組を求めたい。

2変数なので3次元に図示するとこんな感じである。最小値はひとつ定まりそう。真っ赤なところがお盆の低いところになる。

f:id:MikuHatsune:20170502171100p:image

x <- y <- seq(-10, 10, length=300)
f <- function(x, y) x*(x-2)*(x-1)*(x+8) + y*(y-2)*(y-3)*(y+7) -8
xy <- expand.grid(x, y)
z <- f(xy[,1], xy[,2])
i <- 100
cols <- rainbow(i)[cut(z, quantile(z, seq(0, 1, length=i-1)))]
plot3d(xy[,1], xy[,2], z, col=cols, xlab="x", ylab="y", zlab="z")

ちなみにPython の3D プロットは機嫌が悪いようで使えなかった。イラッ☆

離散的に解いた最小値は

xy[which.min(z),]; min(z)
           Var1      Var2
22864 -5.785953 -4.916388
[1] -1245.717

 

その1:解析的に解く

関数が決まっていて、微分が出来て、微分=0 が解析的に解けるならばそれで解いたほうが精確であるし、それが最小値であることが保証される。

高校数学までで言えば、f(x)x関数のとき、x微分して導関数が0となるときのx が得られて、増減表から凹になる部分を探す、というのがセオリー

大学数学でいえば、パラメータがいくつかあるときに、各パラメータ微分偏微分)して偏微分たちの連立方程式が0 となるパラメータセットを考えればよい。

今回のf(x,y)微分できそうな感じで用意したので、微分しよう。微分にはsympy を用いる。

from sympy import *
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
x, y = symbols("x y")
f = x*(x-2)*(x-1)*(x+8) + y*(y-2)*(y-3)*(y+7) -8

parms = [x, y]                           # 変数
deriv = map(lambda z: diff(f, z), parms) # 各変数での偏微分
sol = solve(deriv, parms)                # 偏微分したものを連立方程式で解く

z = map(lambda z: f.subs({x:z[0], y:z[1]}).evalf(), sol)
z = np.array(z, dtype=np.complex)

xy = np.array(sol, dtype=np.complex)

sol はsympy で解いた偏微分方程式オブジェクトで、ところどころに虚数I が入り乱れている。sympy のオブジェクトにはsubs 関数で、辞書形式{x:x0, y:y0} と入れれば、代入したことになる。このとき、sympy 自体がまだ関数であればdoit() をすることで計算され、計算した値が欲しければevalf() で数値を得られる。

sympy 形式でI で表示されている場合、numpy に戻すにはdtype=np.complex を指定する。

すると、

z
array([   11.74548233 +6.77626358e-21j,   -10.78588513 +1.33638236e-51j,
        -565.66379583 +3.34095589e-49j,     4.81250937 +6.77626358e-21j,
         -17.71885808 -8.45930031e-49j,  -572.59676879 -5.13170824e-49j,
        -668.34781483 +6.77626358e-21j,  -690.87918229 +4.37999317e-47j,
       -1245.75709299 +4.41326909e-47j])
xy
array([[ 0.43596395 +2.71050543e-20j,  0.83707381 +1.05879118e-22j],
       [ 0.43596395 +2.71050543e-20j,  2.56096376 -5.29395592e-23j],
       [ 0.43596395 +2.71050543e-20j, -4.89803757 +3.38813179e-21j],
       [ 1.58881707 +5.29395592e-23j,  0.83707381 +1.05879118e-22j],
       [ 1.58881707 +5.29395592e-23j,  2.56096376 -5.29395592e-23j],
       [ 1.58881707 +5.29395592e-23j, -4.89803757 +3.38813179e-21j],
       [-5.77478101 +2.64697796e-23j,  0.83707381 +1.05879118e-22j],
       [-5.77478101 +2.64697796e-23j,  2.56096376 -5.29395592e-23j],
       [-5.77478101 +2.64697796e-23j, -4.89803757 +3.38813179e-21j]])

となり、z をみると最小値-1245 というのがあって、それに対応するxy は(-5.775, -4.898) であり、悪くない。

 

その2:勾配法(急速降下法)

関数の式は立式できても、微分がなかなか解けないという状況はある。そのとき、「変数を少し動かして、もともと最小化したかった関数がより小さくなるようにする」というのを計算機的にゴリ押しする。

アルゴリズム的に書けば、i 番目のパラメータi+1 番目にいい感じに更新するときの具合が、関数から定まる勾配によって支配され、

w^{(i+1)} ¥leftarrow w^{(i)}-¥eta¥nabla E(w^{(i)})

と書ける。ここkで、¥nabla が(偏微分から求まる)勾配、wパラメータセット、¥eta が勾配に応じてどれくらい動かすかの学習率である。

学習率はハイパーパラメータなので、事前にいい感じに設定してやる。アルゴリズムが進化すればこのハイパーパラメータすらも勝手に決まっていく。

¥eta=0.005 とし、(0,0) から始めると、10回くらいの繰り返しで最小値に行っている。この手法は遅い法で、しかも局所解に落ち込むこともある。

解析的に微分できない場合は、微分の定義に戻って、f(x,y) から微小距離¥Delta だけ動かして、¥Delta ¥to ¥infty の極限をとるのを模倣して、¥Delta=10^{-8} くらいで近似すればいい。この方法は、解析的に微分できる場合でも、2通りで試して正解だったらちゃんとできている、というのがわかるので、やっておくとよいらしい。

ここでは、f(x,y)¥to f(x+¥Delta, y+¥Delta)微分ではなく、演算誤差の精度の都合上、本当に見たい点(x,y) を挟んだf(x-¥Delta,y-¥Delta)¥to f(x+¥Delta, y+¥Delta) の差分を見ている。なので割り算するときも2¥Delta である。

# 最小化する関数
def J(x, y):
  return x*(x-2)*(x-1)*(x+8) + y*(y-2)*(y-3)*(y+7) -8

# 勾配
# 偏微分で与えられるならそれでもよい
def gradient_descent(x, y, h=10e-8):
  dx = (J(x+h, y) - J(x-h, y)) / (2*h)
  dy = (J(x, y+h) - J(x, y-h)) / (2*h)
  return np.array([dx, dy])

xy = [0, 0]
np.array(map(lambda z: z.subs({x:xy[0], y:xy[1]}).evalf(), deriv))

def GD(init, eta=0.005, epoch=100):
  res = np.zeros((epoch, len(init)))
  res[0] = init
  cost = np.zeros(epoch)
  cost[0] = J(init[0], init[1])
  for i in range(1, epoch):
    init -= eta * gradient_descent(init[0], init[1])
    res[i] = init
    cost[i] = J(init[0], init[1])
  return res, cost

s1 = GD([0, 0])

plt.subplot(121)
plt.plot(s1[0][:,0], s1[0][:,1])
plt.xlim(-10, 1)
plt.ylim(-10, 1)
plt.title("x and y coordinate")
plt.subplot(122)
plt.plot(s1[1])
plt.title("Value to be minimized")
plt.show()

f:id:MikuHatsune:20170502171136p:image

 

その3:共役勾配法

共役なベクトル列がうまく選べると、急速降下法はその1点の勾配の向きまかせに動いていたが、ちゃんと最適化方向に進行方向が決まるらしい。しかもどれくらい行くべきかの学習率も勝手に決まる。実装力の低いPython novice はライブラリを使う。

scipy にoptimize という一連の最小化最適化ライブラリがあるので、いい感じに使う。書き方はここからパクる。

最小化したい関数J は、パラメータセット¥theta関数にしておくのが簡単である。他に引数が欲しければ、*args と書いておいて関数内で各変数にばらす(後述)。

# conjugate gradient
# 急速降下法と引数の書き方を少し変える
def J(theta):
  x, y = theta
  return x*(x-2)*(x-1)*(x+8) + y*(y-2)*(y-3)*(y+7) -8

def gradient_descent(theta):
  h = 10e-8
  grad = np.zeros(len(theta))
  di = np.diag(np.repeat(h, len(theta)))
  for i in range(len(theta)):
    grad[i] = (J(theta+di[i]) - J(theta-di[i]))/(2*h)
  return grad

init_theta = [10, 10]
res = optimize.fmin_cg(J, init_theta, fprime=gradient_descent, full_output=True)
Optimization terminated successfully.
         Current function value: -1245.757093
         Iterations: 9
         Function evaluations: 29
         Gradient evaluations: 29

この程度の簡単さなら実行速度の恩恵は実感できないが、高速らしい。結果のres には最終的な最小値と、その時のxy の組が入っている。

disp=False にすればコンソールに出力しない、full_output=True にすればすべての出力結果を取得できる。

R ではRcgmin パッケージでできて

library(Rcgmin)
fr <- function(xy){
  x <- xy[1]
  y <- xy[2]
  z <- x*(x-2)*(x-1)*(x+8) + y*(y-2)*(y-3)*(y+7) -8
  return(z)
  }

gr <- function(xy){
  h <- 10e-8
  grad <- rep(0, length(xy))
  for(i in seq(xy)){
    h1 <- replace(rep(0, length(xy)), i, h)
    grad[i] <- (fr(xy+h1) - fr(xy-h1)) / (2*h)
  }
  return(grad)
}

Rcgmin(fn=fr, gr=gr, par=c(10, 10))
$par
[1] -5.774781 -4.898038

$value
[1] -1245.757

$counts
[1] 875 129

$convergence
[1] 1

$message
[1] "Too many function evaluations (> 866) "

 

ようやくここまで来たので、球面調和関数の最適な回転を見つけたい。

球面調和関数の係数は、もとの物体を回転させなくても、係数を回転させればよい。この回転は

c_l^m=¥sum_{m^{’}=-l}^l c_l^{m^{’}}D_{m^{’}m}^l

で得られる。これはsympy ではWigner D 行列としてある。

回転させるパラメータ空間は(¥alpha,¥beta,¥gamma) であり、これらは¥[0,2¥pi¥] の間で連続に変化するはずである。なので、最小値(と最大値)がひとつ定まることは保証されている、はず。ただし局所解であることは多々あると思う。

ここで、最適化関数J には、パラメータ¥theta=(¥alpha,¥beta,¥gamma) を入れ、L^2 ノルムの計算に必要な、球面調和関数係数の各項を*args から入力することにしている。

from sympy import *
from sympy import pi
from sympy.physics.quantum.spin import Rotation
from sympy.functions import re
from scipy import optimize
import numpy as np
import numpy.linalg as LA
import math
from scipy.special import sph_harm as _sph_harm
import pyshtools as p
import matplotlib.pyplot as plt

l, m, n, a, b, g = symbols("l m n a b g") # mp is n
WigD = Rotation.D(l, m, n, a, b, g)
v1 = [10, 1, 2, 0, 1, 5, 3, 4, 0]
r = flatten([[sum(transpose(Matrix(v1[(k**2):((k+1)**2)])) * Matrix([WigD.subs({l:k,m:j,n:n0}) for n0 in range(-k, k+1)])) for j in range(-k, k+1)] for k in range(Lmax+1)])
r = map(lambda z: z.subs({a:abg[0],b:abg[1],g:abg[2]}).doit().evalf(), r)
r = np.array(r, dtype=np.complex) # 完成品

でもこれ、shtools.rotate と結果が合わない…

diff微分もできるけれども、計算がものすごく遅い。

というわけで、数値微分によって勾配を求めたものとする。

def tri(clm_coef):
  vec = clm_coef.coeffs
  return np.concatenate([np.fliplr(vec[1])[:,:-1], vec[0]], axis=1)

def invariantL2vec(vec):
  return np.array([np.linalg.norm(vec[range(k**2, (k+1)**2)]) for k in range(Lmax+1)])

def clm2vec(clm, Lmax=None):
  if Lmax == None:
    Lmax = clm.lmax
  return np.concatenate([tri(clm)[k][(Lmax-k):(Lmax+k+1)] for k in range(Lmax+1)])

lmax = 2
x = p.SHCoeffs.from_zeros(lmax)
ls = [0, 1, 1, 1, 2, 2, 2, 2, 2]
ms = [0, -1, 0, 1, -2, -1, 0, 1, 2]
value = [10, 1, 2, 0, 1, 5, 3, 4, 0]

for (l, m, v) in zip(ls, ms, value):
  x.set_coeffs(v, l, m)

PI = pi.evalf()
a, b, g = PI/2, PI/3, PI/4
x0 = x.rotate(a, b, g, degrees=False)

def J(theta, *args):
  clm0, clm1 = args
  #return LA.norm(clm2vec(clm0) - clm2vec(clm1.rotate(theta[0], theta[1], theta[2], degrees=False)))
  return LA.norm((clm0 - clm1.rotate(theta[0], theta[1], theta[2], degrees=False)).coeffs)

# partial derivatives
def gradient_descent(theta, *args):
  clm1, clm2 = args
  h = 10e-8
  grad = np.zeros(3)
  di = np.diag(np.repeat(h, 3))
  for i in range(3):
    grad[i] = (J(theta+di[i], clm1, clm2) - J(theta-di[i], clm1, clm2))/(2*h)
  return grad

init_theta = np.array([1, 1, 1])
optimize.fmin_cg(J, init_theta, fprime=gradient_descent, args=(x0, x))

Warning: Desired error not necessarily achieved due to precision loss.
         Current function value: 2.529874
         Iterations: 2
         Function evaluations: 17
         Gradient evaluations: 5
array([ 4.49933254, -6.11515311,  3.59954578])

初期値が(1,1,1) だと最小値が0にならなかったようである。いま、球面調和関数の係数x(¥frac{¥pi}{2},¥frac{¥pi}{3},¥frac{¥pi}{4}) だけ回転させたものがx_0 であり、最小値は0となって一致しないとおかしい。

というわけで変な局所解にハマったようである。

本当ならばこういうときに探索する領域をいい感じに飛ばす手法があるっぽいが、適当に初期値を振ることにすると

optimize.fmin_cg(J, np.random.rand(3)*2*pi.evalf(), fprime=gradient_descent, args=(x0, x))
Warning: Desired error not necessarily achieved due to precision loss.
         Current function value: 0.000000
         Iterations: 37
         Function evaluations: 94
         Gradient evaluations: 81
array([1.57079627211901, 1.04719757847306, 0.785398109529357], dtype=object)

となって、確かに(¥frac{¥pi}{2},¥frac{¥pi}{3},¥frac{¥pi}{4}) が推定されている。

ちなみに、x_0 を回してx に一致させるには、(¥alpha,¥beta,¥gamma) で回した時には(-¥gamma, -¥beta,-¥alpha) で逆変換すると戻る。

 

勾配法最適化についてはここらあたり。

パターン認識と機械学習 上

パターン認識と機械学習 上

これは確率的勾配法やモーメント法、AdaGrad がどのような最適化経路をたどるかを可視化しているのでわかりやすい。

 

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証

トラックバック - http://d.hatena.ne.jp/MikuHatsune/20170502/1493712702