본문으로 바로가기



Optimization with R - 1. Newton’s method를 통한 해찾기

category Statistics/Optimization 2017.06.24 22:46
Optimization with R - 1. Newton’s method를 통한 해찾기

\(f(x) = 0\)의 해를 찾아라.

오늘 우리가 공부할 내용은 수치해석학의 가장 유명한(?) 방법인 Newton's method에 대해서 알아보도록 하겠다. Newton's method 방법은 어떤 주어진 함수 \(f(x)\)를 0으로 만드는 입력값(즉, \(f(x)=0\)의 해)을 반복적인 수치계산을 통하여 찾는 방법이다. 이 방법을 적용하기 위해서는 함수 \(f\)가 특정 조건들을 만족해야 하는데, 기본적으로 다음의 두 가지는 꼭 만족해야만 한다.

  • 함수 \(f(x)\)를 0으로 만든다는 것에서, 함수는 결과값이 한 개인 real-value function이어야 한다는 것을 알 수 있다.
  • 우리가 사용하는 함수가 미분가능한 함수, 즉 \(f'(x)\)를 구할 수 있는 함수이어야 한다.

Newton’s method 프로세스

gif from: 위키피디아 Newton’s method

gif from: 위키피디아 Newton’s method

Newton’s method가 어떻게 해를 찾아 가는지는 위키피디아에 링크되어 있는 위의 애니메이션을 보면 도움이 될 것이다. 애니메이션을 보면 현재 우리가 가지고 있는 해의 추정값을 \(x_{old}\)이라고 할 때, 다음 단계의 추정값 \(x_{new}\)를 다음의 단계를 거쳐서 구하는 것을 알 수 있다.

  • step1: \(x_{old}\)를 설정한다.
  • step2: \(x_{old}\)에서의 함수값 \(f\left(x_{old}\right)\)를 구한다.
  • step3: \(x_{old}\)에서의 미분계수 \(f'\left(x_{old}\right)\)를 구한다.
  • step4: 스텝 1~3의 정보들을 가지고, 다음 스텝의 \(x_{new}\)값을 구한다. \[ x_{n+1}:=x_{n}-\frac{f\left(x_{n}\right)}{f'\left(x_{n}\right)} \]
  • step5: 일정 조건을 충족할 때까지, \(x_{new}\)\(x_{old}\)로 바꾼 후 step 1~4를 반복한다.

알고리즘의 직관적인 이해

스텝 4에서 사용한 식을 살펴보면 Newton’s method가 다음 단계의 x를 결정함에 있어서 다음과 같은 몇가지 특징을 가진다는 것을 알 수 있다.

1. 현재 우리가 가진 x가 해일 경우, 움직이지 않는다.

만약 우리가 현재 가지고 있는 \(x_{old}\)가 진짜 우리가 찾고있는 해일 경우, 즉, \(f\left(x_{old}\right) = 0\) 일 경우, 다음의 x값 역시도 그 자리에 머물러 있는 것이 좋을 것이다. 이러한 성질은 아래에서와 같이 newton’s method에 잘 반영되어 있다. \[ x_{n+1}:=x_{n}-\frac{0}{f'\left(x_{n}\right)}=x_{n} \]

2. 현재 함수값이 0에서 멀리 떨어져 있을수록 다음 스텝이 커지고, 현재 자리에서의 기울기가 가파를수록 다음 스텝이 작아진다.

우리가 현재 주어진 \(x_{old}\)에 있다고 생각할 경우, \(x_{old}\)에서의 함수값이 0에서 멀리 떨어져있다는 것은 우리가 찾는 해와 현재 우리가 가지고 있는 \(x_{old}\)와는 멀리 떨어져있을 가능성이 크다. 왜냐하면 현재 0보다 훨씬 크거나 작은 값을 가지는 연속함수를 0으로 만들기 위해서는 현재의 점 \(x_{old}\)을 그만큼 왼쪽으로든 오른쪽으로든 이동해야하기 마련이므로 다음 스텝사이즈는 크게 움직이는 것이 합리적일 것이다.

함수값에 따라서 스텝사이즈를 크게 정했다고 하더라도, 만약 현재 자리의 미분계수의 절대값이 크다면, 그것은 다시 한번 생각해 볼 일이다. 왜냐하면 미분계수의 절대값이 크다는 이야기는 현재의 \(x_{old}\)값을 조금만 움직여도 함수값이 확확 변한다는 의미이기 때문에, 다음 스텝사이즈는 기울기에도 영향을 받아야하는 것이다. 위의 두 가지 속성을 반영하기 위하여 Newton’s method는 \[ \frac{f\left(x_{n}\right)}{f'\left(x_{n}\right)} \] 항을 사용하고 있다.

3. 다음 \(x_{new}\) 방향은 ‘함수값/기울기’ 부호의 반대방향으로 결정한다.

이제 앞에서 살펴본 \(\frac{f\left(x_{n}\right)}{f'\left(x_{n}\right)}\)항을 빼줄지 더해줄지를 결정해야 하는데, 여기서 항을 더 해준다는 것은 현재 위치에서 그 벡터가 가리키는 방향으로 간다는 것이고, 빼준다는 것은 -1을 곱해 더하는 것과 같으므로 현재위치에서 벡터가 가리키는 반대방향으로 간다는 이야기이다.

다음의 필자가 (발로)그린 4가지의 경우를 살펴보자.

함수값과 기울기의 4가지 케이스

그림에서 알 수 있듯, 함수를 0값으로 가깝게 만들려면 \(\frac{f\left(x_{n}\right)}{f'\left(x_{n}\right)}\) 항이 가리키는 방향(검은색 화살표)의 반대방향(빨간색 화살표)으로 점을 움직여야 한다.

Newton’s method의 R로 구현하기

이제까지 우리는 Newton’s method의 식 \[ x_{n+1}:=x_{n}-\frac{f\left(x_{n}\right)}{f'\left(x_{n}\right)} \] 의 의미에 대하여 알아보았다. 이제는 이것을 직접 R로 구현하는 방법에 대하여 알아보도록 하자.

함수는 먼저 해를 찾을 함수 f와 그 함수의 도함수 df, 시작 포인트 x_0, 언제 멈출것인지에 대한 기준값 stop_criteria, 최대 반복횟수 max_iter, 그리고 특정 범위 안에서 찾고 싶은 경우를 대비하여 left_bright_b를 모수로 갖도록 설정하였다.

루프를 스탑하는 기준은 세가지로 설정하였는데,

  • 우리가 설정한 최대 반복횟수가 넘었을때
  • 더이상 x의 값이 움직이지 않는다고 판단되었을 경우
  • 함수 값이 충분히 0에 가까워졌을 경우

세 가지 중 한가지라도 충족이 될 경우 루프에서 빠져나와 결과값을 알려주도록 정의되었다.

newton_method <- function(f, df, x_0,
                          stop_criteria = .Machine$double.eps,
                          max_iter = 10^3,
                          left_b = -Inf, right_b = Inf){

    # criteria for starting point x_0
    if (left_b > x_0 | right_b < x_0){
        warning("your starting point is out of bound,
                please try point in the bounds")
        return(NA)
    } 
    
    # variable initialization
    iters <- 0
    x_1 <- 0
    distance <- Inf
    
    # do while until
    # 1) x does not move
    # 2) exceed the maximun iteration
    # 3) the function value become close to zero
    while ((distance > stop_criteria) &
           (iters < max_iter) &
           (abs(f(x_0)) > stop_criteria)){
        # newton method
        x_1 <- x_0 - ( f(x_0) / df(x_0) )
        
        # out of bound case
        if (x_1 <= left_b) {
            x_1 <- left_b
        }
        if (right_b <= x_1){
            x_1 <- right_b
        }
        distance <- abs(x_0 - x_1)
        iters <- iters + 1
        x_0 <- x_1
    }
    
    # return the final point
    if (x_0 == left_b | x_0 == right_b){
        warning("\n", "Final point is on the boundary", "\n")
    }
    if (iters > max_iter){
        warning("\n", "Cannot find the root during the iteration.", "\n")
    }
    cat("Final value of function:", f(x_0), "\n")
    x_0
}

위에서 적용한 함수가 잘 작동하는지 테스트하기 위해서 cos(x)의 해를 찾아보도록 하자. 먼저 그래프를 그려보도록 하자.

library(ggplot2)
p <- ggplot(data = data.frame(x = 0), mapping = aes(x = x))
p <- p + stat_function(fun= cos, color="blue", size=1.0) + 
         geom_hline(yintercept=0) + xlim(0, 10) + ylim(-1, 1)
p  + ggtitle("Cos(x)")

g <- function(x){ - sin(x)}
newton_method(f = cos, df = g, x_0 = 1, left_b = 0, right_b = 2.5)
## Final value of function: 6.123032e-17
## [1] 1.570796
newton_method(f = cos, df = g, x_0 = 4, left_b = 2.5, right_b = 5)
## Final value of function: -1.83691e-16
## [1] 4.712389
newton_method(f = cos, df = g, x_0 = 6.5, left_b = 5, right_b = 7.5)
## Warning in newton_method(f = cos, df = g, x_0 = 6.5, left_b = 5, right_b = 7.5): 
## Final point is on the boundary
## Final value of function: 0.3466353
## [1] 7.5
newton_method(f = cos, df = g, x_0 = 8, left_b = 7.5, right_b = 10)
## Final value of function: 3.061516e-16
## [1] 7.853982

세번째의 경우는 범위 안에 해가 존재하지 않는 것을 알 수 있다. 다음과 같이 한꺼번에 함수를 적용할 수도 있을 것이다.

roots <- sapply(c(1,4,8), newton_method, f = cos, df = g)
## Final value of function: 6.123032e-17 
## Final value of function: -1.83691e-16 
## Final value of function: 3.061516e-16
roots
## [1] 1.570796 4.712389 7.853982

위에서 구한 해들을 코사인 그래프 위에 표시해 보도록 하자.

# plotting the points on the cosine graph
roots_p <- data.frame(x_star = roots, y_star = cos(out))
p <- p + geom_point(data = roots_p, aes(x = x_star, y = y_star), 
               color = "red", size = 1.5) + theme_bw()
p + ggtitle("Roots of Cos(x) between (0, 10)")

(Special thanks to Sang Woo Ham who help me to puts the roots on the graph, and Seung Hoon Yang who posts the code that inspired me of this post.)

Reference

[1] Wikipedia Newton’s method page: https://en.wikipedia.org/wiki/Newton%27s_method


SHARE TO



티스토리 툴바