6.5 뉴튼-랩슨 알고리즘(Newton-Rhapson Algorithm)

임의의 함수 \(f(x)\)가 주어졌을 때 \(f(x) = 0\) (\(f(x)\)의 해)를 만족하는 \(x\)를 반복적인 수치계산을 통해 찾는 방법

  • Newton-Raphson (N-R) 방법 적용 시 \(f(x)\)의 만족 조건
    • \(x\)의 특정 범위 내에서 \(f(x) = 0\)를 만족하는 유일한 실수값 존재
    • \(f(x)\)는 미분 가능한 함수

Newton-Raphson 알고리즘 예시

N-R 알고리즘(스케치)

  • step 1: 초기치 \(x_{old}\)를 설정

  • step 2: \(x_{old}\)에서 \(f(x_{old})\) 값 계산

  • step 3: \(x_{old}\)에서 접선의 기울기(미분계수) \(f'(x_{old})\) 계산

  • step 4: \(f'(x_{old})\)의 접선이 \(x\)축과 만나는 점을 새로운 값 \(x_{new}\)로 업데이트

\[ x_{new} = x_{old} - \frac{f(x_{old})}{f'(x_{old})} \]

  • step 5: 일정 조건을 만족할 때 까지 step 1 ~ step 4 반복

  • step 4에서 초기값 \(x_0\)이 주어졌을 때 \(f(x_0)\)의 접선은 \(f'(x_0)\) 이고 \((x_0, f(x_0))\)를 통과하므로 접선의 식은 아래와 같음

\[ f(x) = f'(x_0)(x - x_0) + f(x_0) \]

  • \(f(x) = 0\) 일때 \(x\)의 값은

\[ x = x_0 -\frac{f(x_0)}{f'(x_0)} \] - 따라서 다음 단계에서 해의 근사치 \(x_{1} = x_0 - f(x_0)/f'(x_0)\) 이고, 이를 조금 더 일반화 하면,

\[ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} \]

  • 위 식은 테일러 전개(Taylor expansion)를 통해 도출 가능(한 번 생각해 볼 것!!)

N-R 알고리즘의 특징

  1. 현재 \(x_{old}\) 또는 \(x_{n}\)이 0을 만족할 경우, 더 이상 다음 단계로 가지 않음.
  2. 현재 함수값이 0에서 멀리 떨어져 있을수록 다음 스텝이 커지고, 미분계수의 절대값이 클수록 다음 스텝이 작아짐
    • 미분계수의 절대값이 크다 \(\rightarrow\) \(x_n\)을 조금만 움직여도 함수값이 크게 변한다는 의미
    • 따라서 미분계수의 값을 다음 스텝에 반영해야 함.
  3. 다음 \(x_{new}\)의 방향은 \(f(x_{old})/f'(x_{old})\) 부호와 반대방향으로 결정
  4. 수렴속도가 빠르지만 초기값에 따라 알고리즘의 성능이 달라짐
  5. \(f'(x)\)를 반복적으로 계산해야 하고, 경우에 따라 \(f'(x) = 0\)이면 반복식 계산이 불가

반복 종료 조건

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

N-R 알고리즘 구현

  • 알고리즘에 입력되어야 할 변수
    • 초기값과 해를 찾을 범위 지정 \(\rightarrow\) 만약 초기값이 해당 범위를 벗어난 값이 입력되었다면 함수 종료
    • 함수
    • 반복횟수
    • 0과 충분히 가까운 상수(종료 시 필요) \(\rightarrow\) tol
  • 함수 내부 또는 함수 외부에서 1차 미분 함수가 요구
    • 함수 인수로 입력 vs. 함수 내부에서 도함수 계산?
    • 도함수 계산 시 위 예제에서 사용한 R 내장 함수 사용 vs. 미분식 사용?

\[ \lim_{d \rightarrow 0} \frac{f(x + d) - f(x)}{d} \]

  • 반복 종료조건에 도달할 때 까지 반복이 필요 \(\rightarrow\) while 문 사용
  • 반복 조건: 반복이 최대 반복수보다 작고 \(|f(x_{new})|\) 값이 tol 값보다 클 때 까지

newton_raphson <- function(FUN, # 함수
                           x0 = 1, # 초기값
                           max_iters = 5000, # 최대 반복 횟수
                           tol = 1.0e-9, 
                           range = c(-Inf, Inf), 
                           ...) 
{
 iters <- 1;
 grads <- deriv(as.expression(body(FUN)), "x", function.arg = TRUE)
 # grads 반환값 중 "gradient" 값 = f'(x0)
 gap <- x0 - FUN(x0)/attr(grads(x0), "gradient") 
 
 while(iters < max_iters & abs(gap) > tol) {
   # x_new 계산 
   x_new <- x0 - FUN(x0)/attr(grads(x0), "gradient")
   gap <- FUN(x_new) 
   # x_new 가 범위를 벗어난 경우 처리
   if (x_new <= range[1]) x_new <- range[1] 
   if (x_new >= range[2]) x_new <- range[2]
   iters <- iters + 1
   x0 <- x_new # 초기값 업데이트
 }
 
 if (x_new == range[1] | x_new == range[2]) 
   warning("마지막 점이 x 범위의 경계선 상에 있습니다.")
 if (iters > max_iters) 
   warning("최대 반복 때 까지 해를 찾지 못했습니다.")
 cat("x 가", x_new, "일 때 함수값:", FUN(x_new), "\n")
 return(list(solution = x_new, iteration = iters))
 
}

## test: 위 동영상과 동일한 해를 갖는지 확인
f <- function(x) 5 * x^3 - 7 * x^2 - 40 * x + 100 
newton_raphson(FUN = f, 
               x0 = 1,
               range = c(-10, 10)) -> sols
x 가 -3.151719 일 때 함수값: -3.547029e-11