非線形最適化入門(私が)

1.要旨

　今回は，非線形最適化手法の中でも最も有名で，よく利用されているニュートン法について概略をまとめ，RおよびPythonでの実装例を載せました．多変量解析のある分析手法のパラメータ推定で利用した(する予定)ので，勉強も兼ねて実装しました．

2.非線形最適化

2.1 何それ？

　そのまんまですが，非線形関数の最適化(=最小化/最大化)のことです．非線形関数の例としては，$y = x^2$や$y = x\log_e x$などが挙げられます．例えば，$y = x\log_e x$を例にとると，この関数の最小値を与える点を得ようと思った場合，($0<x$で下に凸であることはわかっているとして，)$x$に関して微分して，
$$\frac{d}{dx}y = \log_e x + 1$$
この勾配が$0$となる点の関数の値を計算すれば良いです．つまり，
$$\log x_e + 1 = 0$$
の根を求めれば良いのです．ちなみに$\hat{x} = 1/e$です．
　このように，方程式をスパーンと単純に解ける場合(解析的に解ける場合)は，それでいいのですが，解が陽に表せない場合も多いです．そういった場合に非線形最適化手法が活躍します．

2.2 Newton-Raphson法

非線形最適化の中でも，最も有名で，シンプルで多岐にわたる場面で利用されている手法として，ニュートン法(Newton method, Newton-Raphson method, wiki)があります．数理的な詳細に関しては，こちら(ニュートン法とは何か？？ニュートン法で解く方程式の近似解)をご覧ください．丁寧に図解してくださっていて，高校数学がふんわりわかっていれば，大枠が掴めると思います．これを実装するためには，

*　勾配を計算するための数値微分
*　ニュートン法の更新式の反復

の二つが必要です．一つめに関しては，とりあえずは，

$$\frac{d}{dx}f(x) = \underset{h\rightarrow 0}\lim \frac{f(x + h) - f(x)}{h}$$

を実装しましょう．Rでは，

numerical.differentiate = function(f.name, x, h = 1e-6){
  return((f.name(x+h) - f.name(x))/h)
}

で，pythonでは，

def numerical_diff(func, x, h = 1e-6):
    return (func(x + h)-func(x))/h

みたいな感じでしょうか．これらとニュートン法による逐次計算を用いて，先ほどの

$$\log_e x + 1 = 0$$

を解いてみましょう．Rでの実行例は，

newton.raphson(f.5, 100, alpha = 1e-3) # 0.3678798

pythonでの実行例は，

newton_m(f_1, 100, alpha = 1e-3) # 0.36787980890600075

どちらも，先ほどの方程式の解$1/e\fallingdotseq1/2.718 = 0.3679176$に近くなっていることがわかります．

3.サンプルコード

雑コードですが...

{r, newton_raphson.R}

numerical.differentiate = function(f.name, x, h = 1e-6){
  return((f.name(x+h) - f.name(x))/h)
}

newton.raphson = function(equa, ini.val, alpha = 1, tol = 1e-6){
  x = ini.val
  while(abs(equa(x)) > tol){
    #print(x)
    grad = numerical.differentiate(equa, x)
    x = x - alpha * equa(x)/grad
  }
  return(x)
}

f.5 = function(x)return(log(x) + 1)
newton.raphson(f.5, 100, alpha = 1e-3)

{python, newton_raphson.py}

import math
def numerical_diff(func, x, h = 1e-6):
    return (func(x + h)-func(x))/h

def newton_m(func, ini, alpha = 1, tol = 1e-6):
    x = ini
    while abs(func(x)) > tol:
        grad = numerical_diff(func, x)
        x = x - alpha * func(x)/grad
    return x

def f_1(x):
    return math.log(x) + 1

print(newton_m(f_1, 100, alpha = 1e-3))

ニュートン法が，大域的最適解に収束するには，最適化する関数が，最適化する区間$[a,b]$に大域的最適解が存在し，その区間で凸でなければならないことに注意しましょう．私が昔かいたプログラムは，多変数の場合だったので，理解がいまいちでしたが，1変数の場合だと，シンプルで理解しやすいです．計算の精度は，tolパラメータとhパラメータに依存します．