[ 과제 2 ] Gauss-Newton Method

똑똑한 후배에게 Gradient descent 를 설명해 주었더니, learning rate 구현을 빼먹어서 error function 이 발산하는 쪽으로 가버렸네요..

최적화 방법 중에 그럼, learning rate 없이 할 수 있는 방법을 없을 까요?

*Gauss-Newton Method 설명을 위해서 Newton's Method 에 대한 설명이 필수적이라, 먼저 간략히 설명 뒤에 Gauss-Newton Method에 진행하겠습니다.

Newton's method

Wikipedia : https://en.wikipedia.org/wiki/Newton%27s_method

 - 임의의 시작점에서 함수의 미분을 이용하여 1차 식으로 근사하여 0과 만나는 다음의 x 값으로 x를 업데이트 진행하는 방식입니다. 아래 보시면 gradient descent 와 다르게 learning rate(=step size) 가 보이지 않는 것을 알 수 있습니다.

$f : [ a , b]  \rightarrow \mathbb{R} $    

a,b 범위 안에서 미분가능한 함수 $f(x) = 0$ 근을 반복적인 방법을 통해 찾는 과정 입니다. 

기울기 $f'(x_1)$ 와 한 점 $(x_1,f(x_1))$ 만 알고 있다면 아래와 같이 다음 $x$ 값을 업데이트 할 수 있습니다.

$$ \begin{matrix} f'(x_1) , \quad [ x_1, f(x_1) ] &\rightarrow& f'(x_1)(x-x_1) - f(x_1) = 0 \\ x &=& x_1 - { f(x_1) \over f'(x_1) } \end{matrix} $$

우리가 원하는 오차함수가 운이 좋게 E(x) = 0 이란 값이 나올 수 있지만, 현실적으로 E(x) = 0 이란 값이 나올 수 없는 경우가 더 많습니다. least square method 방법으로  오차함수를 정의 했기 때문에, 잠깐 생각을 바꿔서 미분값이 0이 되는 경우를 찾으면, 이는 분명히 최소, 최대 일 경우 둘중에 한가지로 생각해 볼 수있습니다. 여기에 더해서 오차함수의 특성상 최대값이 나오는 형태가 어려움으로 대부분 최소로 생각을 정리 해볼 수 있습니다. 

다시 한번 설명 드리면 E(x) = 0 의 경우가 존재 하지 않을 수 있기 때문에, E'(x) = 0 으로 문제를 바꿔서 수식을 정리한 것이 아래 내용입니다. 

$$ \begin{matrix} x_{t+1} &=& &x_t&- &{ E'(x_t) \over E''( x_t) }& \\ \rightarrow \boldsymbol{ p_{t+1} } &=& &\boldsymbol{ p_t }& - &H_E(\boldsymbol{ p_t })^{-1} \nabla E(\boldsymbol{ p_t })& \\ \end{matrix} $$

한 단계 더 나아가 matrix 형태로 수식을 정리하면 Gradient, Jacobian 에 이어서 Hessian 이란 무시무시한 놈이 나타납니다. Hessian 의 경우, 2차 미분이라고 간단히 생각을 정리하고 넘어가도록 하겠습니다. Hessian 을 직접 구하지 않고 근사하는 방식이 바로 Gauss-Newton Method 방식이죠!! 

저희가 구현할 내용도 바로 이  Gauss-Newton Method 입니다. 여기까지 왜 Gauss-Newton-Method가 필요한지 설명을 드렸는데,, Hessian 구하기가 그렇게 쉽지 않기 때문이죠!!

Gauss-Newton method

- 접근 방법은 taylor series 를 통해 residual 을 1차 미분 까지만 전개하여 근사화 하여 이를 통해서 파라미터를 업데이트 하겠다는방법이다. 아래 우선 수식을 먼저 공유드리고 세부적인 내용인 차후 업데이트 하도록 하겠습니다. 관심 있으신 분들은 아래 파라미터 업데이트 방식 간단히 python 코드로 구현 해보시기 바랍니다.

 여기서의 핵심은  Hessian 연산을 피하기 위해 테일러 급수 를 통해서 residual 을 근사화하고 이를 통해서 파라미터를 업데이트를 진행하겠다는 방향입니다. 이부분만은 꼭 기억해두시 바래요.

$$ \begin{matrix} \text{Taylor series} \\ f(x) &=&  f(a) + f'(a)(x-a) +{ f''(a) \over 2! }(x-a)^2 \cdots  {  f^{(n)}(a) \over n! }(x-a)^n  &=&  \sum_{n=0}^\infty { f^{(n)}(x-a) \over n!}(x-a)^n \\ \text{approximation} \\ \boldsymbol{r}(\boldsymbol{p}) &\approx& \boldsymbol{r}(\boldsymbol{p_t}) + J_r(\boldsymbol{p_t}) \boldsymbol{r}(\boldsymbol{p - p_t}) \\ \end{matrix} $$

$$ \large \begin{matrix} E(\boldsymbol{p}) = \sum_{i=1}^n { [y_i - f(x_i,\boldsymbol{p})) ] }^2 = \sum_{i=1}^n {r_i}^2 = \boldsymbol{r^T r} \\ \\ \\ { \partial E(\boldsymbol{p}) \over \partial \boldsymbol{p} } = { \partial \boldsymbol{r^T r} \over \partial \boldsymbol{p} } = 2\boldsymbol{r^T} { \partial \boldsymbol{r} \over \partial \boldsymbol{p} } \\ \\ \\ \therefore \boldsymbol{p_{t+1}} = \boldsymbol{p_t} - (J_r^T(\boldsymbol{p_t})J_r(\boldsymbol{p_t}) )^{-1}J_r^T(\boldsymbol{p_t})\boldsymbol{r(p_t)} \end{matrix} $$

**아래 코드를 한번 구현해 보시면 생각이 정리 될 수 있어요

-  그림을 보시면 단 1번의 iteration  으로 아래와 답에 유사하게 나왔습니다. 엄청나죠,,!!

  (기존 naive gradient descent 또는 momentum 에 비교해도 수렴이 확실히 빠르네요)

import numpy as np
import numdifftools as nd
import matplotlib.pyplot as plt

def optimizer_gauss_newton(xdata,ydata,max_iter,epsilon=0.0):

    aidx = 0;
    f, axarr = plt.subplots(2,2)
    axarr[0,0].plot(xdata,ydata,'ro')
    axarr[0,0].set_title('original')

    #model 2th-polynomial: ax^2 + bx +c
    #pt = np.random.rand(3,1)
    pt = np.zeros((3,1))

    r_fun = lambda p:(ydata -(p[0,0]*xdata**2+p[1,0]*xdata+p[2,0]))
    Jrp = nd.Jacobian(r_fun)

    for i in range(1,max_iter):
        r = ydata -(pt[0,0]*xdata**2+pt[1,0]*xdata+pt[2,0])

        pt = pt - np.dot(np.linalg.pinv(Jrp(pt)),r)
        error = np.sum(np.absolute(r))
        if(i%1 == 0 or i == max_iter -1):
            print('iter=[{0}], error={1}'.format(i,round(error,4)))
            aidx +=1
            axarr[int(aidx/2),aidx%2].plot(xdata,ydata,'ro')
            axarr[int(aidx/2),aidx%2].plot(xdata,pt[0,0]*xdata**2+pt[1,0]*xdata+pt[2,0])
            axarr[int(aidx/2),aidx%2].set_title('iter[{0}]'.format(i))
    print(pt)
    plt.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
    plt.setp([a.get_yticklabels() for a in axarr[:, 1]], visible=False)
    plt.show()

if __name__ == '__main__':
    xdata = np.reshape(np.arange(0,1,0.1),(-1,1))
    ydata = 2.1*xdata**2 -1.5*xdata + 0.5
    optimizer_gauss_newton(xdata,ydata,4)