勾配法による最適化 - 驚異のアニヲタ社会復帰の予備

MikuHatsune2017-05-02

ある関数があって、最小化もしくは最大化したい。これらは正負を入れ替えればよいでの最小化を考える。
　
いま、適当に $f(x,y)=x(x-1)(x-2)+y(y-2)(y-3)(y+7)-8$ という関数があったとする。これの最小値とそのときの $(x,y)$ の組を求めたい。
2変数なので3次元に図示するとこんな感じである。最小値はひとつ定まりそう。真っ赤なところがお盆の低いところになる。

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 を指定する。
すると、

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)})$
と書ける。ここｋで、 $\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()

　
その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)$ で逆変換すると戻る。
　
勾配法と最適化についてはここらあたり。